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Abstract 


Constructing  a  model  for  data  in  'R?  is  a  common  problem  in  many  scientific  fields, 
including  pattern  recognition,  computer  vision,  and  applied  mathematics.  Often,  little 
is  known  about  the  process  which  generated  the  data  or  its  statistical  properties.  For 
example,  in  fitting  a  piecewise  linear  model  the  number  of  pieces  as  well  as  the  knot 
locations  may  be  unknown.  Hence  the  method  used  to  build  the  statistical  model 
should  have  few  assumptions  and  yet  still  provide  a  model  that  is  optimal  in  some 
sense.  Such  methods  can  be  designed  through  the  use  of  genetic  algorithms. 

In  this  paper  we  examine  the  use  of  genetic  algorithms  to  fit  piecewise  linear  functions 
to  data  in  R?.  The  number  of  pieces,  the  location  of  the  knots,  and  the  underlying 
distribution  of  the  data  are  assumed  to  be  unknown.  We  discuss  existing  methods 
which  attempt  to  solve  this  problem  and  introduce  a  new  method  which  employs 
genetic  algorithms  to  optimize  the  number  and  location  of  the  linear  pieces.  We  prove 
theoretically  that  our  method  provides  near-optimal  functions  and  present  the  results 
of  extensive  expreiments  which  demonstrate  that  the  proposed  method  provides  better 
results  than  existing  spline  based  methods.  We  conclude  that  our  method  represents 
a  valuable  tool  for  fitting  both  robust  and  non-robust  piecewise  linear  functions. 


Key  Words;  Genetic  Algorithms,  Optimization,  Statistical  Data  Analysis,  Splines, 
Neural  Networks 


1  Introduction 


One  of  the  purposes  of  statistical  data  analysis  is  to  determine  a  functional  relationship 
between  input  and  output  variables  given  a  dataset  of  noisy  observations.  The  dataset, 
denoted  as  D  —  {{xi,yi),  i  =  1,. ,  N},  is  assumed  to  be  a  number  of  realizations  of 
some  underlying  process  combined  with  random  noise,  i.e. 

y  =  f{x)  +  e 


where  e  has  mean  zero.  The  problem  of  determining  f{x)  given  a  set  of  data  points 
and  various  assumptions  is  relevant  to  many  application  fields  including  engineering, 
chemometrics,  and  materials  science  [1,  2]. 

The  construction  of  a  functional  model  for  f{x),  denoted  by  f{x),  can  involve  classi¬ 
cal  statistical  tools  such  as  kernel  methods,  regression,  and  splines  [3,  4],  as  well  as 
more  recent  techniques  such  as  neural  networks,  radial  basis  functions,  and  genetic 
algorithms  [5,  6,  8].  In  using  these  tools  the  function  f{x)  is  often  assumed  to  be  of 
a  certain  form  and  to  meet  certain  mathematical  requirements,  such  as  continuity  or 
differentiability.  In  the  assumed  form  of  f{x)  is  usually  linear,  piecewise-linear,  or 
curvilinear.  Since  the  linear  models  represent  a  subset  of  the  piecewise-linear  models, 
and  any  curve  can  be  approximated  by  a  piecewise-linear  function,  we  chose  to  fit 
piecewise-linear  models  to  our  datasets.  In  this  case,  we  cannot  assume  strict  differen¬ 
tiability.  However,  usually  it  is  the  choice  of  technique,  rather  than  the  data  or  prior 
process  knowledge,  that  motivates  the  placement  of  artificial  mathematical  restrictions 
on  the  final  model  [9].  It  should  also  be  noted  that  not  all  techniques  can  guarantee 
convergence  to  a  near-optimal  f{x).  Given  these  considerations  and  our  desire  to  opti¬ 
mize  the  number  of  pieces  as  well  as  their  placement  in  the  model,  we  utilized  genetic 
algorithms  as  our  primary  model  fitting  tool. 

Genetic  algorithms  (GAs)  are  stochastic  optimization  tools  from  the  field  of  artificial 
intelligence  [11].  They  are  capable  of  finding  near-optimal  solutions  to  problems  with¬ 
out  the  usual  mathematical  model  restrictions  [9,  10].  In  this  work  we  discuss  how 
genetic  algorithms  were  employed  to  fit  piecewise  linear  models  to  data  sets  in  TZ'^, 
where  the  number  of  pieces  (or  knots)  as  well  as  their  placement  are  unknown.  We 
first  present  the  problem  and  discuss  current  methods  for  function  approximation.  We 
then  introduce  genetic  algorithms  and  detail  how  GAs  can  be  used  to  fit  near  optimal 
piecewise-linear  functions.  Several  examples  are  presented  with  results,  comparisons 
are  made  to  existing  methods,  and  areas  for  future  research  are  mentioned. 


2  Problem  Statement  and  Current  Methodology 

Suppose  we  are  given  a  data  set  D  =  {(xi,  yi),  (x2,  y2),  •  ■  ■  ,  (xn,  yN)}  C  TZ^  where  the 
values  (xi,yi)  are  related  by  an  unknown  function  /  such  that  yi  =  /(xj)  -f-Cj,  where  Cj 
is  a  random  error  with  zero  mean  and  constant  variance.  We  assume  that 
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•  /  is  an  /j.*-piecewise  linear  function  where  h*  is  an  unknown  positive  integer 
representing  the  number  of  pieces,  1  <  h*  <  h^a.x,  hmax  known. 

•  The  knot  locations  (zii,zi2), . . .  ,  (z(h*+i)i,  Z(h»+i)2)  of  /  are  unknown  (the  end¬ 
points  of  /  are  also  considered  knots). 

The  problem  is  to  fit  a  function  to  the  data  such  that 

•  The  fitted  function  maximizes  the  quality  of  the  fit  over  all  such  /i-piecewise 
linear  functions,  1  <  h  <  hmax,  where  the  quality  of  the  fit  is  determined  by  a 
specified  function  FF. 

FF  is  defined  in  such  a  way  that  the  fitted  function  is  ‘close’  to  /.  We  make  no 
assumptions  regarding  the  smoothness  of  functions  around  their  knots. 

Classical  approximation  theory  suggests  several  methods  for  fitting  piecewise  linear 
functions,  methods  which  involve  building  models  from  linear  combinations  of  nonlinear 
functions  [12].  Such  linear  estimators  can  be  expressed  as 

f{x)  =  f2Kx{x,Xi)yi  (1) 

where  Kx{x,Xi)  is  a  weighting  function  which  depends  on  some  parameter (s)  A  [3]. 
Some  examples  include  kernel  estimates  and  series  approximators  (which  we  will  not 
explicitly  discuss),  splines,  locally  weighted  regression,  and  the  more  recent  neural 
network  and  radial  basis  function  estimates. 

A  spline  of  order  L  with  knots  at  (zn,  Z12), . . .  ,  (zki,  ZK2)  is  a  function  s  of  the  form 

s(^)  =  H  +  XI  (2) 

2=0  2=1 

for  e  F,  i  =  0, ...  ,L  —  1,  and  5j  G  F,  j  =  1,...  ,K  [13].  In  other  words, 
a  spline  is  a  piecewise  polynomial  where  the  pieces  are  tied  at  knots  in  such  a  way 
that  s  satisfies  certain  continuity  properties  (e.g.,  the  first  L  —  1  derivatives  of  s  are 
continuous).  As  such  they  can  be  viewed  as  an  extension  of  polynomial  regression  [3]. 
Different  classes  of  splines  can  be  formed  by  using  different  basis  functions,  e.g.,  B- 
splines,  periodic  splines,  etc.  Splines  are  useful  when  we  want  an  estimate  which  meets 
a  certain  quality  of  fit  (fitness)  as  well  as  a  smoothness  criterion.  For  example,  we  may 
estimate  y  by  choosing  /  to  minimize 

n-^  '^{y  -  f{xi)f  +  X  f  (3) 

i=i 

for  A  >  0,  m  €  {1, 2, 3, . . . },  and  a  <  Xj  <  6  V  z  =  1, . . .  ,N.  The  solution  /  is  called  a 
smoothing  spline  estimate,  and  A  is  the  smoothing  parameter.  A  determines  the  tradeoff 
between  goodness-of-fit  and  smoothness.  If,  however,  our  fitness  criterion  is  not  of  this 
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form,  a  spline  may  not  be  the  best  type  of  estimate.  Splines  have  applications  in  areas 
such  as  materials  science  [1],  computer  tomography  [3],  and  military  analysis  [14]. 

To  use  (smoothing)  splines  for  analysis,  the  order  L  of  the  spline,  the  number  and 
location  of  the  knots,  the  smoothing  parameter  A,  the  choice  of  basis,  and  the  smooth¬ 
ing  criterion  need  to  be  specified.  This  may  not  be  a  simple  task.  How  to  determine 
‘optimal  values’  for  these  parameters  is  an  area  of  active  research.  Methods  designed 
to  find  a  good  estimate  for  A  are  computationally  demanding  [15]  with  GCV  estimates 
tending  to  oversmooth  the  data  [16].  m  is  often  based  on  prior  information  [3]  as  op¬ 
posed  to  theoretical  considerations.  Schwetlick  and  Schiitze  [19]  describe  an  algorithm 
which  optimizes  the  location  and  number  of  ‘free’  knots  but  is  computationally  ‘too 
expensive’  to  implement  and  involves  the  approximation  of  various  parameters  whose 
effects  on  the  final  estimate  are  unknown.  Larson  [14]  finds  a  closed  form  for  the  min¬ 
imizing  abscissa  for  unknown  knot  locations  of  an  interpolating  spline,  but  does  not 
mention  the  optimization  of  the  number  of  knots.  Suchomski  [20]  and  Kokkonis  [1] 
have  done  some  work  with  fitting  least  squares  variable  knot  interpolating  splines  with 
the  number  of  knots  fixed.  Mammen  [4]  addresses  the  stated  problem  with  locally 
adaptive  regression  splines,  finding  the  best  fit  for  all  A  given  a  fixed  order  L  but  with 
knots  restricted  to  the  class  of  data  points. 

Locally  weighted  regression,  or  LOESS  [21],  fits  an  estimate  g{x)  which  uses  a  neigh¬ 
borhood  of  weighted  observations  whose  values  Xj  are  closest  to  x.  g{x)  may  be  written 
as 

9{^)  =  Y.k{x)yi  (4) 

2=1 

where  theli{x),  i  =  1,. . .  ,N,  depend  on  {xi, . . .  ,  xn},  a  neighborhood  size  q,  a  weight¬ 
ing  function  W,  and  a  distance  measure  p.  Local  fitting  allows  the  estimation  of  a  wide 
class  of  regression  surfaces,  wider  than  that  which  can  be  estimated  by  polynomials. 
However,  LOESS  does  assume  that  y  has  a  normal  distribution  and  that  the  estimate 
g{x)  is  unbiased. 

A  recent  development  in  functional  approximation  is  the  use  of  neural  networks  (NN) 
and  radial  basis  functions  (RBFs).  Multilayer  neural  networks  are  function  approxi¬ 
mators  (in  the  sense  of  (1))  of  the  form,  e.g., 

M 

f{x,W)  =  Y^  Pj9j  (5) 

j=i 

where  Pj,l  <  j  <  M,  are  the  weights  connecting  M  hidden  units  to  the  output  unit, 
oijyl  <j  <  M,  are  weights  connecting  the  input  layer  unit  to  the  jth  hidden  layer  unit, 
and  the  g/s  are  the  hidden  layer  activation  functions  [12].  W  represents  the  matrix  of 
network  weights.  Because  of  the  form  of  /,  splines  and  kernel  estimates  are  sometimes 
viewed  as  special  cases  of  NN  models.  Another  special  case  of  NN  models  is  based  on 
radial  basis  functions  where  the  approximation  is  produced  by  passing  each  Xj  through 
a  set  of  basis  functions  (each  containing  a  RBF  center),  multiplying  the  result  by  a 
coefficient,  and  then  summing  the  results.  In  other  words, 

M 

f{xk,  W)  =  WjT-P(f){\\x  -  CjW/r)  (6) 

i=i 
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where  (f)  is  the  radial  basis  function,  {c^- :  j  =  1, . . .  ,  M}  is  the  set  of  RBF  centers,  r  is 
a  scale  parameter,  and  p  is  the  dimension  of  the  data.  Often  (j)  corresponds  to  a  Gaus¬ 
sian  density  [5,  6].  Note  that  a  radial  basis  function  network  (RBFN)  is  essentially  a 
kernel  method  for  regression.  In  general,  NNs  are  easily  programmed  and,  as  a  result, 
have  become  an  almost  universal  optimization  ‘crank’:  simply  toss  in  the  data,  add 
any  number  of  parameters,  and  wait  for  gradient  descent  to  produce  the  result  [22]. 
However,  there  is  no  method  for  finding  the  best  network  architecture  for  a  given  prob¬ 
lem.  Once  an  architecture  has  been  chosen  (heuristically),  the  system  often  requires 
numerous  adjustments,  supplied  by  an  experienced  user,  and  do  have  a  tendency  to 
overfit  or  overparameterize  the  data  [22].  They  may  also  get  stuck  at  local  minima, 
unlike  GAs  [5,  25].  Ghen  and  Jain  [12]  report  that  backward  propagation  can  be  slow 
and  sensitive  to  noise.  They  suggest  a  robust  modification  whose  parameters  are  the 
focus  of  further  study.  RBFNs  have  been  shown  to  outperform  multilayer  perceptrons 
(MLPs)  [6]  even  though  the  choice  of  centers  and  the  curse  of  dimensionality  can  make 
implementation  difficult.  It  should  also  be  noted  that  both  NN  and  RBFN  results  lack 
interpret  ability  [22]. 

A  discussion  of  the  advantages  and  disadvantages  of  the  above  methods  led  us  to 
consider  the  possibility  of  using  genetic  algorithms  (GAs)  as  a  data  fitting  tool.  GAs 
need  three  things  to  construct  a  model  for  a  dataset:  a  way  to  calculate  a  response,  an 
error  function  (i.e.,  a  measure  of  quality  of  fit),  and  a  set  of  parameter  values.  However, 
unlike  splines,  the  values  given  to  the  parameters  are  not  critical  to  the  algorithm’s 
result  or  its  convergence  -  they  only  influence  the  rate  at  which  the  algorithm  converges. 
By  utilizing  GAs  we  may  discard  the  artificial  requirements  of  differentiability  and 
smoothness  imposed  on  the  functional  form  of  the  model  by  the  above  methods  [9].  It 
is  also  possible  to  define  the  fitness  (or  optimization)  function  in  such  a  way  that  the 
same  algorithm  can  be  used  to  fit  both  robust  and  non-robust  functions  -  robust  in  the 
sense  of  being  insensitive  to  outliers.  In  the  case  of  piecewise  linear  models,  the  power  of 
GA  optimization  allows  not  only  variable  knot  placement  but  also  a  variable  number  of 
knots.  Such  variable  models  have  been  fit  using  Bayesian  estimation  as  opposed  to  GA 
optimization  [23,  24],  but  GAs  can  guarantee  convergence  to  an  optimal  solution  for  a 
sufficiently  large  number  of  iterations  [25].  Related  works  include  those  of  Karr  [2,  8], 
who  has  applied  GAs  to  the  problem  of  least  squares  (LS)  and  least  median  squares 
(LMS)  curve  fitting  where  the  specific  form  of  the  curve  is  known  (e.g.,  a  polynomial 
of  degree  2)  and  the  data  is  noiseless,  and  Vankeerberghen  [26],  who  have  used  GAs  to 
obtain  the  parameters  for  specific,  known  laboratory  system  models  (e.g.,  a  hyperbolic 
model)  which  are  nonlinear  in  the  parameters.  We  intend  to  design  a  GA  based  method 
which  fits  both  robust  and  non-robust  near-optimal  piecewise  linear  functions  to  data 
in7^2. 


3  Genetic  Algorithms 


Genetic  algorithms  are  stochastic  search  methods  which  provide  a  near  optimal  so¬ 
lution  to  the  evaluation  function  of  an  optimization  problem  [9].  They  can  be  used 
to  search  complex,  multimodal  surfaces  via  steps  which  have  been  designed  to  mimic 
the  processes  of  natural  genetic  systems.  They  work  simultaneously  on  a  number  of 
possible  solutions  to  help  prevent  the  algorithm  from  getting  stuck  in  a  local  optimum. 
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With  little  adjustment  -  in  most  cases  by  only  redefining  the  fitness  function  -  GAs  can 
be  applied  to  a  wide  range  of  optimization  problems.  Thus  it  is  possible  to  use  the  same 
basic  algorithm  to,  for  example,  fit  lines  satisfying  different  optimization  criterions  to 
the  same  dataset.  Situations  in  which  the  effectiveness  of  GAs  has  been  demonstrated 
can  be  found  in  the  literature  of  various  scientific  fields  such  as  scheduling,  classifier 
systems,  and  pattern  recognition  [10]. 

In  GA  algorithms,  each  possible  solution  is  encoded  as  a  string  or  chromosome;  a  set 
of  such  chromosomes  is  called  a  population.  An  evaluation  (fitness)  function  provides 
a  mapping  from  the  chromosome  space  to  the  solution  space.  The  algorithm  starts 
with  an  initial  population  of  a  fixed  number  of  randomly  generated  strings.  At  each 
iteration,  three  basic  operations  -  selection,  crossover,  and  mutation  -  are  applied  over 
the  current  population  to  yield  a  new  population  of  strings.  This  cycle  is  repeated  until 
some  termination  criterion  is  met,  at  which  time  the  best  string  achieved  is  generally 
taken  as  the  solution  to  the  optimization  problem. 

We  will  briefly  outline  the  basic  framework  of  a  GA  before  discussing  specific  GA 
concepts  that  can  be  applied  to  our  problem. 


3.1  Basic  Framework 

For  a  more  detailed  look  at  the  GA  framework,  we  will  discuss  the  stages  of  a  specific 
GA  model,  the  elitist  model.  Consider  the  problem  of  maximizing  a  function  FF(x)  G 
V,  where  V  is  a  finite  set  and  FF(x)  >0^  x  eV.  Each  string  S,  built  from  members 
of  a  finite  alphabet  A  =  {cki,.  . .  ,o;a},  corresponds  to  a  value  x  in  V  and  may  be 
written  as 

S  —  ('Yi,  '72)  •  •  •  )  ^igth)i  '^i  ^  A  Vi  1,...  ,  Igth 

The  number  of  different  strings  that  are  possible  is  For  our  work  we  will  be  using 
binary  strings,  i.e.,  A  =  {0, 1},  so  a*®*'*  =  2*®*'*.  A  random  sample  of  size  M  (M  even) 
is  drawn  from  these  possible  strings  with  replacement  to  form  the  initial  population 
Q.  The  evaluation  or  fitness  value  of  each  string  S  is  fit(S)  =  FF(x)  where  x  G  F  is 
the  value  represented  by  S. 

The  first  operation,  selection,  is  modeled  after  Darwin’s  concept  of  ‘survival  of  the 
fittest’.  Strings  from  the  population  are  selected  and  placed  in  a  mating  pool;  the 
probability  of  selection  for  string  j  is  =  fit(Sj)lY^^\fit{Si).  For  example,  if 
PO)  =  X)i_i  M  strings  are  selected  and  placed  in  the  mating  pool  by  repeating 
the  following  process  M  times: 

1.  Generate  a  random  number  rndi  from  [0,1] 

2.  If  rndi  <  select  5i;  for  j  =  2, . . .  ,  M,  if  <  rndi  <  P^,  select  Sj 

Note  that  strings  with  low  fitness  values  are  rarely  selected  while  some  strings  may  be 
selected  more  than  once.  The  mating  pool  is  taken  as  our  new  population  Qi. 

In  single  point  crossover,  or  reproduction,  pairs  of  strings  exchange  information,  thereby 
generating  two  new  offspring  for  the  next  population.  All  strings  are  paired  at  random 
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in  such  a  way  that  each  string  belongs  to  only  one  pair  (hence  there  are  M/2  pairs). 
If  the  given  pair  is  denoted  as 

P  =  (/?i, . . .  ,  pigth)  and  r  =  (n, . . .  ,  Tigth) 

and  Pc  is  the  probability  that  a  given  pair  of  strings  undergoes  crossover,  then  the 
crossover  operation  may  be  described  as 

1.  Generate  a  random  number  rnd  from  [0,1] 

2.  If  rnd  >  Pc,  no  crossover  is  performed.  If  rnd  <  Pc, 

(a)  Generate  a  random  integer  pos  from  [l,lgth-l]. 

(b)  Strings  ^  and  r  are  replaced  by  strings  jd'  and  r'  where 

=  (^1,  .  •  •  ,  Pposi  '^pos+\i  •  ■  •  1  Tigth)  and  T  =  (Ti,  .  .  .  ,  Tpps,  Ppos+lj  ■  ■  ■  i  Plgth) 

Once  this  operation  has  been  applied  to  each  pair  of  strings,  the  resulting  population 
is  denoted  Qa-  Q2  has  M  strings,  some  of  which  may  have  also  been  elements  of  Qi. 

Mutation  involves  the  random  altering  of  characters  in  the  strings  of  Q2.  Let  pm  denote 
the  probability  of  mutation  of  a  given  character.  Then,  for  each  character  A  of  every 
string  /3,  the  mutation  stage  consists  of 

1.  Generate  a  random  number  rnd  from  [0,1] 

2.  If  rnd  >  Pm,  Pi  is  not  mutated.  If  rnd  <  Pm, 

(a)  If  Pi=0,  then  Pi  becomes  1;  else,  Pi  becomes  0. 


The  mutation  probability  may  vary  over  iterations,  initially  taking  a  high  value,  then 
decreasing  to  a  pre-specified  minimum,  then  increasing  again  in  the  later  stages  of  the 
algorithm  (see  [29]).  When  the  algorithm  has  little  knowledge  of  the  search  space,  it  is 
encouraged  to  explore  its  domain  through  a  high  mutation  probability.  As  the  number 
of  iterations  increases  the  algorithm  will  move  towards  a  solution,  hence  the  mutation 
probability  is  decreased  to  allow  a  search  in  the  vicinity  of  this  solution.  To  avoid  the 
attraction  of  the  algorithm  to  a  local  optimum,  the  mutation  probability  is  increased 
in  the  later  stages  to  again  allow  for  a  broader  search.  The  resulting  population  we 
denote  as  Q3. 

Note  that  through  mutation,  a  given  string  can  become  any  of  the  2*^'*'^  possible  strings. 

We  now  replace  our  initial  Q  with  Q3  and  repeat  the  above  stages  until  the  algorithm 
converges  to  a  satisfactory  solution.  The  stages  we  have  discussed  so  far  are  common 
to  all  GA  models.  In  the  elitist  model  of  GA’s  (EGA),  a  further  operation,  elitism, 
is  added  to  ensure  that  as  the  algorithm  is  running  knowledge  about  the  best  string 
obtained  so  far  is  preserved.  In  this  way  the  algorithm  can  report  at  any  time  the  best 
solution  which  has  been  achieved.  The  basic  steps  which  are  added  to  the  above  model 
to  form  the  elitist  model  are 


7 


1.  Once  the  initial  population  Q  has  been  generated,  find  the  fitness  value  of  each 
string  S  in  Q.  Let  the  string  with  the  maximum  fitness  value  fitmaxQ  of  all  of 
the  strings  in  Q  be  denoted  as  SmaxQ- 

2.  Compare  the  fitness  value  of  each  string  in  Qs  with  the  fitness  value  of  SmaxQ-  If 
no  string  in  Qz  has  a  fitness  value  greater  than  or  equal  to  fitmaxQ,  replace  the 
worst  string  in  Qz  with  SmaxQ- 

At  this  point,  Qz  is  taken  as  the  new  Q  and  the  algorithm  begins  another  iteration  by 
returning  to  the  first  operation,  selection. 

This  basic  framework  can  be  modified  to  specifically  address  the  problem  of  interest. 
One  such  modified  algorithm  which  includes  concepts  which  can  be  applied  to  our 
problem  is  the  variable  length  genetic  algorithm  (VLGA).  These  concepts  will  be  used 
to  design  an  algorithm  for  our  problem  whose  computational  complexity  is  less  than 
that  of  the  VLGA. 


3.2  Variable  Length  Genetic  Algorithm  (VLGA) 

The  variable  length  genetic  algorithm  was  designed  by  Bandyopadhyay,  Murthy,  and 
Pal  [27]  for  optimizing  the  number  of  hyperplanes  needed  to  classify  patterns  in  a 
multidimensional  feature  space.  It  is  able  to  consider  solutions  with  varying  numbers 
of  planes  by  using  strings  of  varying  lengths.  For  example,  a  set  of  h  hyperplanes, 
1  <  h  <  hmax,  hmax  known,  could  be  represented  by  the  string 

S  =  (7i,  72,  •  ■  •  ,  7/i*l);  V?  =  l,...  ,h*L 

where  A  is  the  set  of  possible  characters  and  L  characters  (e.g.,  bits)  are  used  to 
represent  each  plane.  Hence  the  string  length  is  not  fixed  but  variable  -  for  a  given 
h  e  H  =  {1, . . .  ,  hmax},  the  string  length  is  h*  L.  The  algorithm  included  the  basic 
stages  -  selection,  crossover,  and  mutation  -  as  well  as  elitism,  modified  to  handle 
strings  of  different  lengths.  The  elitism  stage  ensures  that  the  algorithm  will  converge 
to  an  optimal  solution,  theoretically. 

The  fitness  criterion  was  defined  so  that  the  set  of  hyperplanes  with  the  maximum 
fitness  value  would  (1)  have  the  minimum  number  of  misclassifications  of  any  set  of  h 
hyperplanes,  heH,  and  (2)  use  as  few  hyperplanes  as  possible  in  doing  so.  Hence  the 
fitness  function  could  be  represented  as 


fit{S)  =  {N-  misss)  +  (7) 

^max 

where  N  is  the  number  of  data  points,  misss  is  the  number  of  points  misclassified  by 
the  set  of  hyperplanes  represented  by  S,  hs  is  the  number  of  hyperplanes  represented 
by  5,  and  hmax  is  the  maximum  possible  number  of  hyperplanes.  Hence  maximizing 
the  fitness  function  will  yield  a  set  of  hyperplanes  which  is  parsimonious  yet  minimizes 
the  number  of  misclassifications. 
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In  designing  the  proposed  algorithm  we  included  an  elitist  stage  ,  as  was  done  in  the 
VLGA,  to  ensure  algorithm  convergence  to  an  optimal  solution.  The  more  important 
concept  borrowed  from  the  VLGA,  however,  is  the  form  of  the  fitness  function.  The 
VLGA  fitness  function  is  based  on  two  criteria  -  minimization  of  the  number  of  mis- 
classifications  and  the  number  of  hyperplanes.  Similarly,  we  are  interested  in  finding  a 
solution  which  (1)  is  ‘closest’  to  the  majority  of  the  data  points  yet  (2)  does  so  using 
the  least  number  of  lines  possible.  Hence  the  form  of  the  proposed  fitness  function  will 
mimic  that  of  the  VLGA  fitness  function. 


3.3  Convergence 

Bhandari,  Murthy,  and  Pal  [25]  have  proven  theoretically  that  an  elitist  genetic  al¬ 
gorithm  (fixed  length  strings)  will  converge  to  an  optimal  string  as  the  number  of 
iterations  goes  to  infinity.  The  two  characteristics  that  are  necessary  and  sufficient  for 
algorithm  convergence  are 

•  The  optimal  string  from  the  present  population  has  a  fitness  value  no  less  than 
the  fitness  values  of  the  optimal  strings  from  the  previous  populations. 

•  Each  string  has  a  positive  probability  of  going  to  an  optimal  string  within  any 
given  iteration. 

By  including  the  elitist  stage  in  the  proposed  algorithm,  we  will  be  able  to  ensure 
convergence  to  an  optimal  solution.  This  convergence  is  independent  of  the  choice 
of  values  for  the  algorithm  parameters  (M,  Pc,  Pm)  although  these  parameter  values 
do  influence  the  rate  of  convergence.  There  is  no  theory  to  indicate  the  number  of 
iterations  necessary  for  convergence.  Two  popular  heuristic  stopping  rules  are 

•  Execute  the  process  for  a  fixed  number  of  iterations  and  report  the  best  string 
found  as  the  solution. 

•  Execute  the  process  until  the  fitness  value  does  not  show  adequate  improvement 
over  a  fixed  number  of  iterations,  and  report  the  best  string  found  as  the  solution. 

We  will  employ  the  first  stopping  rule  in  our  GA  based  method. 


3.4  Piecewise  Linear  Fitting 

We  have  discussed  the  elitist  GA  and  its  convergence  to  an  optimal  string.  This 
convergence  led  us  to  include  elitism  in  our  GA  based  method  for  finding  a  near- 
optimal  solution  to  the  problem  of  fitting  an  h-piecewise  linear  function  to  a  given 
dataset  D  'va.T^.  Designing  such  an  algorithm  required  us  to  ‘encode’  each  possible 
solution  as  a  string,  define  the  fitness  function  to  reflect  the  idea  of  a  ‘properly’  fitted 
function,  and  select  values  for  the  algorithm  parameters  {M,Pc,  etc.).  How  these  steps 
were  achieved  will  be  discussed  in  Sections  4  and  5.1.1. 
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We  have  yet  to  define  the  meaning  of  ‘a  set  of  lines  representing  a  distribution’  -  this 
will  be  done  in  Section  4.  This  meaning  is  critical  -  if  the  given  datapoints  are  from 
this  distribution,  then  the  representative  set  of  lines  is  the  ‘ideal’  solution  which  we 
seek.  Note  that  the  definition  of  any  fitness  function  will  be  dependent  upon  the  given 
dataset  D.  Thus,  for  the  same  dataset  size  N,  we  will  get  different  optimal  solutions 
depending  upon  the  choice  of  D  (assuming  that  all  the  steps  mentioned  above  have 
been  achieved)  regardless  of  the  specific  form  of  the  fitness  function.  However,  if  the 
datasets  are  from  the  same  distribution,  the  results  obtained  by  our  method  should 
be  related.  The  nature  of  the  relationship  between  these  results  needs  to  be  explored 
in  terms  of  the  fitness  function,  theoretically.  Lastly,  the  fitness  function  needs  to  be 
defined  in  a  way  that  as  the  number  of  iterations  and  the  datset  size  N  go  to  infinity, 
the  solution  obtained  will  go  to  the  ‘ideal’  solution. 

Sometimes  the  given  dataset  may  contain  ‘outliers’  which  we  do  not  want  to  influence 
our  determination  of  the  fitness  of  a  given  string.  Note  that  the  meaning  of  ‘outlier’  is 
unclear,  yet  we  would  like  to  define  the  proposed  fitness  function  in  such  a  way  that  the 
‘outliers’  do  not  contribute  to  its  value.  We  shall,  for  the  sake  of  convenience,  assume 
that  at  most  crit  percent  of  the  data  points  are  ‘outliers’  and  remove  them  from  the 
fitness  calculation.  The  value  1  -  crit  is  referred  to  as  the  critical  level.  The  value 
of  crit  depends  on  the  given  dataset  and/or  the  user.  If  the  dataset  has  no  ‘outliers’, 
then  crit  may  be  taken  as  0. 

The  GA  based  optimization  method  to  be  proposed  can  be  used  for  any  given  value  of 
crit.  Hence  this  optimization  scheme  may  be  described  as  robust  from  the  point  of 
view  of  not  yielding  to  the  influence  of  outliers  [26,  28].  Theoretically,  the  properties 
of  the  fitness  function  -  where  the  fitness  function  is  defined  with  the  intention  of 
removing  crit  outliers  -  are  to  be  studied. 

The  present  article  deals  with  all  of  these  theoretical  issues.  In  the  next  section,  the 
phrase  ‘a  set  of  lines  representing  a  continuous  distribution’  is  defined  mathematically. 
A  function  FF,  which  will  later  be  taken  as  our  fitness  function,  is  also  defined.  We 
will  be  using  the  polar  representation  of  a  straight  line  for  the  purpose  of  coding.  This 
polar  representation  will  also  be  used  in  the  development  of  our  mathematical  theory. 


4  Theory  of  Piecewise  Linear  Fitting  in 

4.1  Mathematical  Formulation 

Our  stated  goal  is  to  fit  piecewise  linear  functions  to  datasets  in  Ji?.  In  order  to  do 
this,  we  need  to  mathematically  specify  exactly  what  it  means  for  a  set  of  lines  to  ‘fit’ 
a  datset.  Generally  speaking,  if  we  think  of  a  given  dataset  as  a  realization  of  some 
random  variable,  we  would  like  our  set  of  lines  to  represent  the  center  of  the  density  of 
that  random  variable.  The  statement  of  this  idea  in  mathematical  terms  is  as  follows. 

Let  us  assume  that  we  have  h*  lines  where  1  <  h*  <  fimax)  hmax  is  a  known  positive 
integer.  For  each  i,  i  =  1, . . .  ,h*,  let  the  equation  of  the  i*^  line  be 
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X  cos  6oi  +  y  sin  doi  =  doi 


for  some  Oqi  G  (0,  tt]  and  d^i  G  Tl.  Assume  that  the  and  (i  +  1)*'^  lines  intersect  at 
the  point  (z(i^.i)i, Z(i+i)2).  We  denote  the  first  knot  as  (zii,zi2)  and  the  last  knot  as 

(Z(h*+l)l)  Z(h»+1)2)- 

Let  Co  >  0  and  let  the  set  Bi  be  expressed  as 


fj  =  |(a;,2/)  : 


y  G 


doi  -  X  cos  6oi  doi  -  x  cos  Ooi 

- ^0) - IT-;; - r  Co 


sm( 


sinOoi 


]  X  G  [zji,  Z(i^x)i] 


We  denote  Uf=i  Bi  as  B.  The  probability  density  function  on  B  is  denoted  by  a  :  B 
[0,  oo),  where 


a 


ai{x,y) 

0 


X  G  [zii,z(i+i)i],  i  =  ,h* 

otherwise 


We  define  our  probability  measure  P  as  P{A)  =  a{x,  y)  dx  dy  for  all  Borel  A  C 
B{B),  the  Borel  a-field  of  B. 

Let  (Xi,  Yi),  (X2,  Y2),  ■  •  •  ,(Xn,Yn)  be  independent,  identically  distributed  random 
vectors  with  density  a,  i.e.,  there  exists  a  probability  space  (O,  A,  Q)  such  that  (Xi,  Yi)  : 
(O,  A,  Q)  {B,  B{B),  P),  i  =  ,N,  where  Q(C)  =  P((Xi,  Yi)(C'))  V  C  G  A 

We  also  assume  that 


1.  For  each  i  =  1,. . .  ,h*, 


^  ^  [Zil  > 


(i.e.,  ai  is  symmetric  about  the  line  xcos^oi  +  ysin^oi  =  doi) 

2.  ai{x,  y)  >  0  y  (x,  y)  e  Bi  y  i  =  1, . . .  ,  h* 

3.  a(x,  y)  and  ai{x^  y),  i  =  1,...  ,  h*,  are  continuous. 

4.  Jb  a(x,  y)  dxdy  =  l 


(8) 


In  the  above  mathematical  formulation,  the  primary  assumption  is  the  assumption  of 
symmetry  of  the  underlying  density  around  a  piecewise  linear  function  (assumption 
#1).  The  other  stated  assumptions  are  ordinary  properties  of  a  continuous  probability 
density  function. 

We  will  fit  a  piecewise  linear  function  to  a  given  dataset  by  utilizing  a  GA  to  search 
for  the  function  described  above.  Hence  we  need  to  represent  our  solution  space  in  a 
way  compatable  with  the  search  method  of  a  GA. 
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4.2  Solution  Space 


We  are  given  a  data  set  D  =  {(xi,yi),  (x2,y2),  •  •  •  ,(xN,yN)}  which  is  a  realization 
of  the  random  vectors  (Xi,Yi),  i.e.,  Xi(a;)  =  Xj,  Yi(a;)  =  ji  for  some  w  G  i  = 
I,...  ,N.  Define  X(i)  =  min{xi;  i  =  ,N},  X(jv)  =  max{xi;  i  =  ,N},  y(i)  = 

min{yi;  i  =  1,...  ,N},  y(Ar)  =  max{yi;  i  =  1,...  ,N}.  GAs  try  to  find  an  optimal 
solution  over  a  finite  solution  space.  Thus  the  solution  space,  i.e.,  the  collection  of  all 
h-piecewise  linear  functions  where  h  £%,  H  =  {1, . . .  ,  hmax},  needs  to  be  discretized. 
To  formulate  the  problem  mathematically,  we  proceed  in  the  following  way. 

Let  Lj  denote  the  straight  line 


X  cos  dj  +  y  sin  9j  =  dj 


where  j  belongs  to  an  index  set.  Depending  upon  the  given  context,  the  notation  Lj 
may  also  refer  to  the  function 


Lj{x) 


dj  —X  cos  9 j 
sin  9j 


Let  Lj^.^h  represent  an  h-piecewise  linear  function  with  pieces  Lj,2,h,  ■  ■  ■  ,Lj^h,h) 

where  denotes  the  straight  line  of  the  set  of  h  straight  lines,  i  =  1,...  ,h,  and 
each  Lj^i^h  satisfies  the  following  properties: 

1.  Ljj^h  represents  the  straight  line  x  cos  6 =  dj^i^h  where  9jj^h  (0  < 

<  tt)  is  the  polar  angle  formed  when  the  polar  axis  is  taken  as  the  y-axis  and 
the  origin  is  taken  as  the  intersection  point  between  the  y-axis  and  dj^i^h 

is  the  perpendicular  distance  of  the  line  from  the  point  (0,0). 

2.  For  every  i,  i  =  1,...  ,h  -  1,  and  intersect  and  the  point  of 

intersection  is 

3.  =  X(i),  Zj,(/i+i),/i,i  =  X(jv),  <  2j^(i+i),h,i  V  z  =  1, . . .  ,  h 

4.  Lj^.ji  =  Lj  i^fi  if  ^  ^  ^  ^  1, . . .  fh. 

The  knot  locations  ((zj,i,h,i,  Zj,i,h,2),  •  •  •  ,  (zj,(h+i),h,i,  Zj,(h+i),h,2))  of  any  1  <  h  < 

hmax,  are  not  restricted  to  the  datapoints  {(xi,yi),  (x2,y2), .  ■ .  ,  (xN,yN)}- 

For  each  h,  1  <  h  <  hmax,  let  Ch  represent  the  class  of  all  h-piecewise  linear  functions 
Lj^.^h  which  satisfy  the  above  properties.  Then  Uh^n  is  the  collection  of  functions 
under  consideration. 

Note  that  {jheu^h  is  uncountable.  In  order  to  obtain  a  ‘near  optimal’  function 
from  this  space,  we  need  to  ‘suitably’  discretize  \JheH^hi  i-e-,  larger  discretizations 
should  lead  to  finer  representations  of  the  solution  space.  We  can  achieve  this  by 
restricting  the  values  of  9  and  d.  Let  L  be  the  number  of  bits  used  to  express  9 
and  let  la  be  the  number  of  bits  used  to  express  d.  We  restrict  9j^i^h  to  the  values 


12 


{^>  •  •  •  )  specifying  dj^i^k,  we  utilize  the  rectangle  rect  formed  by 

the  points  (x(i), (a;(jv), 2/(1)),  (a;(i), 2/(Jv))  and  (a;(iv),2/(iv))-  Note  that  rect  contains 
the  entire  given  data  set.  Let  diag  be  the  length  of  the  diagonal  of  rect  and  let  £0  be 
defined  as 


/)  _  (  X(i)  cos  0  +  y(i)  sin  0  ifO<0<7r/2 

^  ~  I  X(jv) cos 9  +  y(i) sin 6  if 7r/2  <9<'k 

Then  for  a  given  dj^i^h  may  only  take  values  within  the  set  {dj^i^h  =  : 

k'j,i,h  G  {0, 1, . . .  ,  2’*^  —  1},  d  =  diag  I  —  1)}.  Note  that  a  line  with  d  =  le  intersects 
rect  at  the  point  (a;(i),?/(i)),  if  0  <  0  <  7r/2,  or  the  point  (a;(jv),  y(i)),  if  7r/2  <  9  <  tt 
(the  value  kj^i^hS,  0  <  kj^i^hd  <  diag,  is  sometimes  referred  to  as  the  offset  value). 

For  each  h,  1  <  h  <  hmax  let  JC)  denote  the  finite  set  of  functions  in  Ch  which 

satisfy  these  restrictions,  where  9  has  ©  possible  values  and  k  has  K  possible  values 
{note  that  0  =  2'*^  and  K.  =  2*'*).  Then  C\{Q,K)  may  be  expressed  as 

Cl{Q,  K)  =  e  Ch  :  Lj^i^h  is  of  the  form 

X cos 9j^i^h  +  y sin 9j^i^h  =  +  hi,h^  i  =  l,...  ,h, 

0j,i,h  e  . . .  ,  tt},  kj,i,n  e  {0, 1, . . .  ,  2'<‘  -  1},  and  5  =  diag/ [2^^  - 1)}. 

Hence  {£°(0,/C)  :  U  =  1,2,...}  and  {Cl{Q,JC)  :  Id  =  1,2,...}  represent  increas¬ 
ing  sequences  of  nested  sets.  For  sake  of  clarity,  we  will  henceforth  specify  Lj^.^h  ns 
L{ej^.,h,kj,.,h)  and  Lj^.^h{x)  as  L{0j^.^h,kj,.,h){x).  Note  that  dj^.^h  and  kj^.^h  represent  the 
vectors  {9j,i,h,  •  •  •  ,  9j,h,h)  and  (kj^i^h,  ■  ,  kjxh)- 

For  L{ej,.,h,  kj^.^h)  e  >C°(0,  /C)  we  define 

E,,a,b{L{0j,.,h,  kj^.^h))  =  {(a;,  y):ye  {L{0j,,h,  kj,.,h){x)  -  e,  L{0j^.,h,  kj,.,h){x)  +  e),  x  e  [a,  6]} 
for  a  <C.  b,  e  >  0,  and  let  Ef^{L{0j^.^fi,kj^.jff)  =  fej^.^/i)). 


If  we  recall  the  piecewise  linear  function  which  lies  at  the  center  of  the  density  of 
{x,y)  and  we  denote  it  as  /,  i.e.. 


d{\i  X  cos  Ooi 
sin  $oi 


0 


X  e  [zii,  Z(j+i)i],  for  i  =  1, . . .  ,  h* 
otherwise 


(9) 


then  the  support  B  of  the  density  is  equivalent  to  £'£o,zii,Z(ft.+i)i(/)- 
Corresponding  to  E^{L{0 j^.^htkj^.^h))  we  define  the  set  of  piecewise  linear  functions 
£<«>  =  €  C„  :  %,»)))  >  0.95} 


where  P{A)  =  f^a(x,y)dxdy  for  A  Borel,  A  c  B{B).  C]^’  represents  those  functions 
L{0j^.^h,kj,.,h)  G  £/i  for  which  the  probability  of  the  region}  (a;,  y)  :  y  G  {L{0j^.^h,,kj^.^h){x)- 
e,L{0j^.^h,kj^.^h){x)  +e),  x  e  [zn,Z{h*+i)i]}  is  greater  than  0.95.  Just  as  {£^(0,/C)  : 
la  =  1, 2, ... }  and  {£^(0,  /C)  :  la  =  1, 2, . . . }  represent  increasing  sequences  of  sets,  we 

have  4"'  c  £)«>  if  <  €2.  These  definitions  involving  E^{L{0j^.^h,  kj^.^h))  will  be  used 
to  define  our  optimization  criterion. 
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Figure  1;  Support  B  =  Bi{]B2  with  center  lines;  Example  function  from  £2(^5 

Figure  1  shows  an  example  dataset  where  the  support  B  —  Bi\J B2  oi the  density  a  is 
centered  around  a  2-piece  linear  function.  A  function  from  £3  also  shown. 

This  completes  an  outline  of  the  mathematical  context  of  our  problem.  We  now  look 
at  the  function  to  be  optimized. 


4.3  Optimization  Criterion 


The  function  which  is  chosen  to  represent  a  given  dataset  is  often  that  which  minimizes 
the  sum  of  the  squared  error  (i.e.,  the  least  squares  function).  However,  a  criterion 
based  on  least  squares  often  yields  a  solution  that  is  not  robust  -  outliers  in  the 
dataset  can  pull  the  least  squares  function  away  from  most  of  the  data  points  and 
hence  away  from  /  [26,  28].  For  this  reason  we  will  not  to  use  least  squares  as  our 
optimization  criterion.  Instead,  we  note  that  our  concept  of  a  ‘fitted’  function  is  one 
which  represents  the  center  of  a  symmetric  density  function.  Given  a  dataset  D,  if 
represents  our  ‘fitted’  function  then  the  majority  of  the  data  points  in 
D  should  fall  within  the  region  {{x,y)  :y  e  {L{d  ikj,.,h>){x)  - + 
e),  X  e.  for  some  ‘small’  e  >  0.  This  observation  is  the  basis  of 

our  optimization  criterion. 


Let  V’j,-,/!  denote  an  h- piecewise  linear  function  where  each  piece  * 

satisfies  the  properties  listed  for  For  any  e  >  0,  define 


1 

0 


l<e 

otherwise 


1,...  ,/i. 


The  proposed  function  to  be  optimized  (i.e.,  fitness  function)  may  then  be  stated  as 


N 


h 


hr. 


■) 


(10) 
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Note  the  similarity  in  form  to  the  fitness  function  of  the  VLGA  (Eqn.  7).  Our  solution 
space  £°(0, /C)  =  some  0  and  JC  represents  the  collection  of  func¬ 

tions  under  consideration.  By  the  continuity  of  a{x,  y),  with  =  U/jgw  we  know 
there  exists  some  e*  :  ^  0  and  =  0  V  e  <  e*  (i.e.,  there  exists  some  function 

L{6j^.^fi,  kj^.^h)  ^  for  some  value  he  71  for  which  P{E^{L{9j^.^h,  kj,.,h)))  >  0.95  when 
e  >  e*  but  no  such  function  exists  when  e  <  e*).  To  balance  accuracy  and  robustness, 
we  specify  that  at  least  95%  of  the  data  points  in  D  (not  necessarily  100%)  fall  within 
e*  of  the  optimal  function  -  our  solution  belongs  to  This  is  equivalent  to  setting 
crit  =  5%  (see  Section  3.4).  The  choice  of  95%  is  heuristic  -  it  can  be  altered  de¬ 
pending  upon  the  characteristics  of  the  dataset  (e.g.,  the  percentage  of  outliers)  and 
the  desired  precision.  Hence  our  goal  is  to  use  genetic  algorithms  to  find  an  optimal 
h'-piecewise  linear  function  kj*,.,h')  ^  -^'’(0,  K)  fl  such  that 


FFe*,N{L{0j*,,h',  kj*,,h'))  =  ^  (11) 


where  h'  —  min{/i  G  Ti  :  ^  ^  0}- 

We  observe  that  if  we  require  that  the  entire  dataset  D  fall  within  e*  of  the  optimal 
function  then,  for  0  and  K,  sufficiently  large,  as  the  number  of  iterations  and  N  ^  oo 
the  optimal  function  should  correspond  to  the  least  squares  function.  This  will  be 
discussed  in  more  detail  in  the  following  section. 

For  the  sake  of  clarity,  a  function  will  be  referred  to  as  optimal  if  (1)  (1  -  crit)%  of 
the  data  points  fall  within  e*  of  the  function  and  (2)  of  all  functions  for  which  this  is 
true,  the  function  contains  the  least  number  of  lines.  Note  that  this  is  different  from 
an  optimal  solution.  An  optimal  solution  is  the  optimal  string  to  which  an  elitist  GA 
converges  and  it  represents  the  string  which  maximizes  the  given  fitness  function.  A 
function  will  be  referred  to  as  properly  fitted  (or  proper)  if  at  least  (l-crzt)%  of  the 
data  points  fall  within  e  of  the  function  for  a  specified  e  >  0  (hence  optimal  functions 

are  also  proper).  Note  that  represents  the  set  of  proper  functions  in  £^.  Unless 
otherwise  specified,  crit  will  be  taken  as  0.05. 

Our  goal  can  be  attained  if  and  only  if 


1.  Our  search  space  £”(0,  K.)  contains  an  optimal  function  as  0  — oo  and  K.  oo. 

2.  An  optimal  solution  represents  an  optimal  function,  i.e.,  a  function  such  that  at 
least  95%  of  the  data  points  fall  within  e*  of  the  function  and  of  all  functions  for 
which  this  is  true,  the  function  contains  the  least  number  of  lines. 

3.  The  proposed  algorithm  converges  to  an  optimal  solution. 


Statement  #3  will  be  discussed  in  Section  4.5.2.  Statements  #1  and  #2  will  be  justified 
in  detail  for  the  case  hm&x  =  1  and  these  results  will  be  extended  to  the  case  of 
hmax  >  1,  h*  unknown,  h*  eTL. 
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4.4  Case  h^ax  =  1 

4.4.1  Optimal  Function  in  Search  Space 

In  the  case  where  hmax  =  1,  h*  and  h'  are  both  known  to  be  1. 

Recall 

£?(0,  /C)  =  i,i)  G  Cl  :  L{0j,i^i,  is  of  the  form 

a;  cos  1,1  +  y  sin0j,i,i  =  +  kj^i^iS,  where 

0j,i,i  is  one  of  0  values,  kj^i^i  is  one  of  K.  values} 

Let 

Vi  =  {L{9m,i,i,km,i,i)  e  £i  :  3  {xr,yr)  e  rect  satisfying  L(0^,i,i,  A:^,i,i), 

0  <  0m, 1,1  ^  ^  T^}- 

Vi  represents  the  set  of  all  lines  in  Ci  which  intersect  rect. 

We  will  justify  the  following  statements: 

1.  For  any  fixed  e  >  e*  our  class  £j(0,/C)  will  contain  a  proper  line  (i.e.,  a  line  in 
C^f^)  as  0  — >•  oo  and  JC  oo. 

11.  For  any  fixed  e  >  e*  the  set  of  proper  lines  in  £°(0,/C)  increases  to  the  set  of 
proper  lines  in  "Pi  as  0  — >  oo  and  K  -)■  oo. 

Additionally,  we  will  show 

III.  For  0  <  p  <  1  let  .4(^)(p)  =  {L{0m,i,i,  km,i,i)  G  Vi  :  P{E,{L{0m,i,u  km,i,i)))  >  p}, 
and  let  e*  :  A^^p\p)  ^  0  and  A^^\p)  =  0  V  e  <  e*.  Then  for  p  =  1.0  and 
N  oo,  A^^p^ip)  converges  to  a  unique  optimal  function  L(^j«,i,i,  /i:j*,i,i)  G  Pi 

As  the  discretization  of  Ci  becomes  finer  (i.e.,  as  0  — oo  and  K,  oo),  Statement  I 
says  that  for  any  e  >  e*  our  solution  space  will  contain  a  proper  function,  hence  when 
e  =  e*  our  solution  space  will  contain  an  optimal  function.  Similarly,  Statement  II 
indicates  that  for  any  e  >  e*  our  solution  space  will  contain  all  proper  functions  (hence 
all  optimal  functions)  which  intersect  rect.  Statement  III  says  that  as  our  critical  level 
increases  to  1.0  (i.e.,  as  we  require  that  a  larger  percentage  p  of  the  data  points  fall 
within  e*  of  the  ‘fitted’  function)  and  N  ^  oo,  the  set  of  optimal  functions  intersecting 
rect  decreases  to  a  single  function.  This  unique  function  is  the  least  squares  function. 

Note  that  we  have  restricted  ourselves  to  considering  only  those  functions  which  in¬ 
tersect  rect  so  our  solution  space  is  actually  £5(0,/C)  restricted  to  those  lines  which 

intersect  rect  (and  similarly  for  Ci^).  This  is  validated  by  the  following  theorem: 
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Theorem  4.1  For  any  e  >  e*  and  for  0  and  K  sufl&ciently  large,  there  exists  an  proper  line 
e  Cl  such  that 

{{x,y)  :  y  =  x  €  [zn,Z2i]}f]  rect  ^  ^ 


A  proof  of  Theorem  4.1  and  a  graphical  representation  of  the  proof  are  provided  in  the 
appendix. 

Prom  Theorem  4.1  we  know  that  if  a  proper  function  exists,  then  there  exists  a  proper 
function  which  intersects  rect.  It  follows  that  without  loss  of  generality,  given  e  >  e* 
we  can  restrict  ourselves  to  ©  and  K.  sufficiently  large  and  only  those  proper  lines  which 
intersect  rect,  i.e.,  those  proper  lines  which  belong  to  Vi- 

We  now  justify  Statements  I,  II,  and  III  with  the  help  of  the  following  propositions 
and  theorems,  the  proofs  of  which  can  be  found  in  the  appendix. 

For  simplicity,  let  p  =  ig  +  k6  for  given  6  and  k. 

Statement  I  will  be  justified  if  we  can  show  that  for  ©  and  JC  large  and  any  e  >  e*, 
given  any  proper  line  in  Vi  we  can  find  a  line  in  £?(©,  /C)  which  is  arbitrarily  close  to 
it.  This  is  implied  by  Proposition  4.1  which  states  that  given  any  line  in  Vi  we  can 
find  a  ©  and  JC  such  that  there  exists  a  line  in  £?(©, /C)  which  is  arbitrarily  close  to 
the  given  line. 


Proposition  4.1  Let  L(0m, i,u  km, i,i)  €  Pi.  Let  ^  >  0.  Then  3  (©^, /C^)  such  that 
V  ©  >  ©^  and  K  >  JC^,  3  L{9,  k)  for  which 

1.  L{e,k)  E  C\{Q,K) 
and 

2.  \  e-  Om,!,!  |<  ^/2  and  \p-  pm, i,i  \<  C/2. 


Now,  given  C  >  0,  we  would  like  there  to  exist  a  ©^  and  JC^  such  that  for  any  ©  > 
©^,  JC  >  JC^,  given  any  proper  line  in  Pi,  £?(©, /C)  will  contain  a  line  which  is 
arbitrarily  close  to  the  given  proper  line.  We  know  such  a  {Q^,JC^)  exists  by  the 
following  theorem. 

Theorem  4.2  For  each  C  >  0,  3  (©^,  JC^)  :  for  all  ©  >  ©^  and  for  all  JC  >  JC^, 
given  any  L{9m,i,i,km,i,i)  €  Pi,  3  L{9,k)  : 

1.  L{9,Ji)  E  £?(©,/C) 
and 
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2.  \  e-  Omxi  l<  ^/2  and  \p-  pmx^  |<  ^/2. 


Thus  given  0  and  K,  sufficiently  large  and  any  e  >  e*  we  can  get  a  line  which  is 
arbitrarily  close  to  a  proper  (or  optimal)  line.  Hence  Statement  I  is  justified. 

For  the  purpose  of  justifying  Statement  II,  let  {0i  :  0j  <  0i+i,  i  =  1,2,...}  and 
{Ki  :  Ki  <  ^i+i,  i  =  1, 2, . . . }  represent  the  possible  values  of  0  and  K.  For  any  e  >  0 
we  define 


and 


and  let  Aei  =  Ai(0.95)  and  A  =  A (0.95).  Ai  represents  the  set  of  all  proper  lines 
which  belong  to  £i(0i,  /C*)  while  A  represents  the  set  of  all  proper  lines  which  belong 
to  Vi.  We  would  like  Ai  ->  as  i  — oo.  However,  we  first  need  that  if  0  ^  oo 
and  K,  oo  then  any  proper  line  which  is  the  limit  of  any  sequence  of  functions  in 
V\  is  contained  in  Vi  (and  hence  belong  to  our  search  space).  This  is,  in  fact,  true,  as 
stated  in  Proposition  4.2. 


Proposition  4.2  For  each  z  =  1,2, . . .  ,  let  T(^ni,i,i)  ^  /^i(0i,A^i)  : 

^  ^lim  nnd  Alnj,l,l  ^  Mim  SOme  ^lim)  0  <  ^Hm  ^  i  nnd 
0  <  him  <  oo,  as  ?  — >•  oo.  Let 

71  =  {L(6>iini,  fciim)  :  3  a  sequence  {T(0„.,i,i,  A:„.,i,i)}^i,  7/(6»„.,i,i,  A:„,,i,i)  e  £5(0i,/Cj) 

such  that  and  km, i,i  him} 

Then  given  any  e  >  e*  the  proper  lines  in  Vi  are  the  proper  lines  in  Ti- 


This  implies  that  Vi  contains  any  proper  (or  optimal)  function  which  is  a  limit  of  a 
sequence  of  functions  in  Vi-  We  may  now  justify  Statement  II. 

Note  that  the  fitness  function  :  (0,7r]  x  [-M,  M]  -)■  [0,  oo)  is  continuous  where 

for  any  line  in  7},  0  belongs  to  (0,  tt]  and  the  distance  of  the  line  from  the  origin  is  less 
than  M.  With  continuous  and  bounded,  we  state  Theorem  4.3. 

Theorem  4.3  Let  Ad,  i  =  1, 2, . . .  ,  and  A  be  as  defined  above.  Then 
Ad  ^  .Ae  as  i  oo. 


Hence  the  proper  (or  optimal)  lines  which  intersect  rect  are  in  the  search  space  as 
0  — >•  oo  and  K  oo. 

The  above  theorem  holds  for  any  Aeip),  0.95  <  p  <  1.0,  and  for  any  ^  ^  {^p}p=0.95) 
where  e*  :  A;(p)  'A  0  and  A(p)  =  0  V  €  <  e*  (i.e.,  if  we  require  that  96%  of  the  data 
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points  in  a  given  dataset  fall  within  e  of  our  fitted  function,  then  ^5.96  is  the  smallest 
e  >  0  for  which  such  a  function  exists).  With  this  in  mind,  we  conclude  with  a  theorem 
which  justifies  Statement  III. 

Theorem  4.4  Let  Ae{p)  and  e*  be  as  defined,  0.95  <  p  <  1.0.  Then  .4,£»^(1.0) 
converges  to  a  single  optimal  line  as  N  ^  oo. 

A  graphical  representation  of  this  theorem  is  provided  in  the  Appendix. 

Suppose  we  require  that  all  of  the  data  points  fall  within  ej  q  of  the  ‘fitted’  function. 
Then  for  0  and  )C  large,  a.s  N  ^  oo  our  optimal  function  will  converge  to  the  least 
squares  function.  This  implies  that  the  proposed  algorithm  can  be  used  to  fit  non- 
robust  as  well  as  more  robust  optimal  functions. 

By  justifying  the  above  statements  we  have  shown  that  our  search  space  contains  an 
optimal  function.  It  remains  to  be  shown  that  the  optimal  solution  of  the  proposed 
GA  is  an  optimal  function  for  /imax  =  1- 


4.4.2  Optimal  Solution 


For  any  e  >  e*  we  know  from  the  above  section  that  we  may  choose  0  and  K  sufficently 
large  so  that  a  properly  fitted  function  is  contained  in  the  search  space 


As  defined  previously,  let 


_  /  1  \y-  L(0j,i,x, 

0  otherwise 


y))  I 

Since  h  takes  only  one  value,  we  may  drop  the  term  (1 


kj,i,i)ix)  \<  e 

T^)  from  Eqn.  10  and  allow 

/imax  ^  ^ 


We  write  as  opposed  to  FFg^N  to  emphasize  the  dependence  of  the  optimal 

solution  on  0  and  K. 


Define 


Note  that 

liniAT-^oo  AFF^^N,e,KiL{9j^i^i,  —  PiUb^JaKb  Fe,a,biL{9j,l,l,  kj^i^i)))  —  P{E^{L[9j^i^i,kj^i^-\))) 


Given  the  convergence  of  the  elitist  GA,  we  know  that  by  including  elitism  in  the 
proposed  algorithm  we  can  ensure  that  it  will  find  an  optimal  solution  L{9j^i^i,  kj^i^i) 
for  a  given  jOP{Q,)C),  e,  N,  and  a  sufficiently  large  number  of  iterations.  So  let 
be  such  that 
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i)6£»(0,x:)  ^-P’e,Ar,0,x:(-C'(^j,l,l>  ^j,l,l)) 


where  the  dependence  of  an  optimal  solution  on  e  and  N  has  been  made  explicit. 

also  maximizes  AFFe,Ar,0,A:(.£'(%,i,i)  To  determine  whether 

the  algorithm  to  be  proposed  converges  to  an  optimal  function,  we  examine 
limAr_+oo^-P'T;,iv,0,x:(T(0j»,i,i,%M,i)e,Ar)  for  e  =  e2jv  where 

A:j*,i,i)£2jv,Jv)  >  0-95.  If 

is  at  least  0.95,  then  the  proposed  algorithm  does  indeed  converge  to  an  optimal  func¬ 
tion  (i.e.,  the  optimal  solution  represents  an  optimal  function).  Theorem  4.5  states 
that  this  is  indeed  true. 


Theorem  4.5  Let  e2N  be  as  defined  above.  Then  for  appropriate  0  and  JC 
liminfAr_^oo^T’Fe2jv,jv(T(0j.,i,i,^j.,i,i)e2^,iv)  >  0.95 


The  proof  of  Theorem  4.5  is  presented  in  the  Appendix. 

Hence  we  have  established  that  for  h^ax  =  !> 

1.  Our  search  space  £*’(0,  /C)  contains  an  optimal  function  as  0  — >  oo  and  IC  co. 

2.  The  optimal  solution  represents  an  optimal  function,  i.e.,  a  function  such  that  at 
least  95%  of  the  data  points  fall  within  e*  of  the  function  and  of  all  functions  for 
which  this  is  true,  the  function  contains  the  least  number  of  lines. 


4.4.3  Remarks 

1.  For  large  0ig  and  /Cig,  |  9i—6i-i  |  and  |  k—ki-i  \  are  both  small,  so  an  optimal  line 

in  £°(0io,/Cjo)  will  be  close  to  an  optimal  line  L{9m*,i,i,km*,i,i) 
in  Vi  in  terms  of  Euclidean  distance. 

2.  A  search  procedure  based  on  the  above  mathematical  formulation  may  be  de¬ 
signed  so  that  0,  JC,  and  e  are  adjusted  adaptively,  i.e.,  in  a  way  that  is  depen¬ 
dent  upon  the  algorithm’s  result.  For  example,  we  may  start  with  initial  choices 
010,  /Cjp,  and  eig,  and  run  the  algorithm  for  a  finite  number  of  iterations.  If  we 
reach  several  proper  results,  we  may  reduce  eig\  if  we  reach  a  single  proper  result, 
we  may  increase  0  and  1C.  We  then  repeat  this  process  until  a  ‘satisfactory’ 
solution  has  been  achieved,  i.e.,  a  ‘fitted’  function  which  reflects  the  center  of  the 
probability  density  function  of  the  observed  random  variables  with  an  acceptable 
level  of  precision.  Note  that  if  0  and  (or)  K  are  (is)  small  or  e  is  large,  it  is  pos¬ 
sible  for  the  algorithm’s  result  to  be  close  to  an  optimum  in  terms  of  probability 
but  not  in  terms  of  Euclidean  distance.  Since  appropriate  values  for  0,  /C,  and  e 
are  unknown  apriori,  we  need  to  implement  an  adaptive  procedure. 
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3.  Recall  that  our  requirement  that  95%  of  the  data  points  be  within  e*  of  the 
optimal  solution  is  heuristic.  Depending  upon  the  particular  dataset  at  hand 
and  the  desired  accuracy  of  the  solution,  we  may  alter  this  criterion  to  better 
suit  our  needs.  For  example,  we  may  choose  our  critical  level  to  be  97%  if  our 
dataset  has  no  outliers  or  90%  if  our  dataset  has  several  outliers. 

4.  In  Statement  III  and  its  corresponding  proof  (Theorem  4.4)  we  justified  that  with 

a  critical  level  of  1.0  and  e  =  e*  the  algorithm  converges  to  a  unique  optimum 
as  — >■  oo  for  a  sufficiently  large  number  of  iterations.  Note  that  this  unique 

optimum  corresponds  to  the  least  squares  solution  (which  is  nothing  but  the  line 
represented  by  Eqn.  9  for  h*  =  1).  Hence  our  method  can  be  used  to  fit  both 
robust  (using  an  e  criterion)  and  non-robust  (using  a  least  squares  criterion) 
solutions. 

5.  It  is  possible  for  our  optimization  problem  to  have  more  than  one  optimal  solu¬ 
tion.  For  example,  consider  Figure  2  which  shows  a  cross-section  of  the  density 
of  the  observed  random  vector  (X,  Y).  Suppose  we  have  two  parallel  lines  lying 
in  the  x  —  y  plane,  one  located  at  the  center  of  the  interval  marked  with  a  number 
1  (crossing  the  a;-axis  at  the  star  marked  with  a  number  1)  and  one  located  at 
the  center  of  the  interval  marked  with  a  number  2  (crossing  the  x-axis  at  the 
star  marked  with  a  number  2).  The  percent  values  on  the  graph  indicate  the 
percentage  of  the  data  points  with  an  x  value  lying  within  the  corresponding 
marked  interval.  Using  these  percentages  we  see  that  95%  of  the  data  points  fall 
within  e  of  line  #1  and  95%  of  the  data  points  fall  within  e  of  line  #2.  Hence 
both  of  these  lines  would  satisfy  our  optimization  criteria  (Eqn.  11). 


Figure  2:  The  existence  of  two  optimal  lines  for  a  given  data  set 


Z 


6.  We  have  assumed  that  the  support  of  ai{x,y),  B,  is  rectangular  in  shape.  How¬ 
ever,  B  may  have  curved  symmetric  boundaries  as  opposed  to  straight  lines.  For 
example,  let  f{x)  be  as  defined  in  Eqn.  9,  i.e.. 


fix)  = 


dm—x  cos  0ni 
sin  $01 

0 


X  e  [2:11,2:21] 
otherwise 
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Then  we  could  have  B  as  shown  in  Figure  3.  All  of  the  above  results  except 
Theorem  4.4  hold  for  such  a  support  B  as  long  as  the  symmetricity  of  ai{x,y) 
is  maintained.  The  symmetricity  is  essential  due  to  the  nature  of  For 

Theorem  4.4  to  hold  we  also  require  that  71  7^  72  (so  that  the  center  portion  of 
/  has  positive  width). 


Figure  3:  A  dataset  whose  density  has  a  support  B  with  curved  boundaries 


4.5  Case  /imax  >  1 

Having  justified  Statements  #1  and  #2  for  the  case  /imax  =  1,  we  now  consider  the  case 
of  >  1,  h*  unknown,  and  h*  G  %.  We  initially  extend  the  results  of  Section  4.4 
to  the  case  of  hmax  >  I,  h*  known,  and  h*  gH. 


4.5.1  Optimal  Function  in  Search  Space 


Suppose  h=  h*,  h*  >  1,  h*  gH  where  h*  is  known. 


Recall  that  our  search  space  is 


Cl*  (0,  K,)  =  e  Ch*  :  Lj^i,h*  is  of  the  form 

X cos +  y sin V  i  =  ,h*, 

Oj,iM  hi,h*  e{0,l,...  ,  5  =  diag/ {2^^ -1)}. 

where  each  L{6j^i^h*,kj^i^h*),  i  =  1,  ■ . .  ,h*,  satisfies  the  properties  listed  in  Section  4.2 
for  Lj,i,h. 


With  respect  to  optimization  we  can  approach  each  ^■s  we  did  the  /imax  =  1  case 

with  ai{x,y)  as  ai{x,y)  and  [zn,  Z(i+i)i]  as  [2:11,^12].  We  then  recognize  £^.(0, /C)  as 
simply  Dill  Uj  £^,ij(0,  a:)  where 
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•  is  £i(0,/C)  restricted  to  [zj,i^h\uZj,i+i,h\i] 
and 

•  the  union  is  taken  over  all  possible  knot  locations  and  over  all  possible  pieces  (any 
knot/piece  combinations  which  contain  pieces  that  do  not  intersect  are  removed). 

By  the  results  of  Section  4.4.1,  for  any  i,  j  we  can  get  arbitrarily  close  to  any  line  over 
[zj,i,h*,i,  Zj,i+i,h\i]  intersecting  rect,  hence  we  can  get  arbitrarily  close  to  any  piecewise 
function  with  knot  locations  whose  x  coordinates  are  [zj^i^h*,i,  ■  ■  ■  ,Zj,{h*+i),h\i]  and 
whose  pieces  intersect  rect.  By  taking  the  union  over  all  possible  knot  locations  with  x 
coordinates  \zj^i^h*,\i  •  ■  ■  i  and  over  all  we  see  that  an  optimal  h*-piecewise 

function  is  contained  in  (0,  /C)  as  0  — oo  and  K,  — >  oo. 

By  defining  our  search  space  as  £°(0,/C)  =  £°(0)  (see  Section  4.3)  we  can 

guarantee,  using  the  same  reasoning  as  in  Section  4.4.1,  that  an  optimal  /i*-piecewise 
function,  h*  unknown,  h*  G  %,  will  be  contained  in  the  search  space  as  0  ^  oo  and 
Ki  — y  oo. 

To  show  that  the  optimal  solution  is  a  properly  fitted  function,  we  first  consider  the 
case  of  h*  known. 


4.5.2  Optimal  Solution 


If  h*  is  known,  h*  G  then  for  appropriate  0  and  /C  the  proposed  algorithm  is 
searching  for  a  function  G  £^.(0, /C)  ^  which  maximizes 

N 

i=l 


)(xi,yi)  +  (l- ^)  (12) 

'^max 
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over  all  L{0j^.^h*ikj,.,h*)  €  i2^.(0,/C).  As  was  done  in  Section  4.4.2,  we  may  drop  the 
term  (1  —  jr^).  The  convergence  of  the  elitist  GA  combined  with  the  continuity  of 
FF^^n  ensures  that  Theorem  4.5  holds,  i.e., 

\immiN-^ooAFF^^^^NiL{0j*,;h*,kj*,.,h*)e2N,N)  >  0-95 


Hence  the  optimal  solution  does  represent  a  properly  fitted  function.  Note  that  for 
any  fixed  h,  the  algorithm  will  find  an  optimal  function  in  \  Hence  an 

algorithm  to  find  an  optimal  function  in  the  /imax  >  1,  h*  unknown,  h*  e'H  case  could 
be  designed  as  follows: 


1.  Divide  >C°(©,  /C)  into  its  component  classes  >C5(0,  /C),  £2(®)  ^))  •  •  ■  >  ^Vax(®> 


2.  On  each  class  £°(0,  /C),  h  =  1, . . .  ,  hmax,  use  an  elitist  GA  to  find  kj*,.,h) 

where  FF^*^^^Q^K{L{0j*,.,h,  FFe*,N,0,KiL(0j,;h,  kj,.,h))- 

3.  Define  Lopt  =  . . .  ,  T(0j*,,w,  Then  the  solution  is 

taken  as  the  function  L{0j*^.^h',kj*,.,h')  ^  Topt  where 


FF,,,N,0,KWj*,,h>,kj*,,h'))  =  max  FF,*,N,0,K{F{Oj*,;h,kj*,,h)) 


This  type  of  GA  we  call  a  partitioned  genetic  algorithm.  We  have  noted  previously 
that,  theoretically,  the  elitist  GA  will  converge  to  an  optimal  solution  for  a  sufficiently 
large  number  of  iterations.  It  is  easy  to  show  (although  the  proof  will  not  be  given 
here)  that  the  proof  of  convergence  to  an  optimal  string  holds  for  this  algorithm  as  well. 
This  justifies  Statement  #3  (see  Section  4.3);  having  previously  justified  Statements 
#1  and  #2,  we  conclude  that  our  research  goal  can  be  achieved.  Due  to  efficiency 
considerations,  we  plan  to  implement  a  partitioned  GA  instead  of  a  VLGA.  However, 
the  optimal  string  of  the  partitioned  GA  will  represent  an  optimal  function  -  by  the 
same  reasoning  as  above  and  by  its  relationship  to  the  variable  length  genetic  algorithm. 

Using  a  partitioned  GA,  for  each  fixed  value  of  /i  G  'H  we  find  the  string  which  maxi¬ 
mizes 

N 

i.e.,  the  minimal  number  of  points  have  distance  greater  than  e*  from  the  function 
represented  by  the  chosen  string.  From  the  resulting  group  of  strings  we  then  select 
as  our  optimal  string  the  one  which  represents  the  least  number  of  lines.  We  observe 
that  this  optimal  string  will  represent  the  function  L{0j*^.^h',  kj*,.,h'  )  G  C^{Q,}C)  which 
maximizes 


N 


FFe*,N,0,K{HGj,;h'^kj^. 


,h'))  = 

i=l 


h' 


hr] 


Noting  the  similarity  between  the  fitness  function  of  the  partitioned  GA  and  the  fitness 
function  of  the  VLGA  (see  Section  3.2,  Eqn.  7),  we  conclude  that  an  optimal  string 
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from  the  partitioned  GA  will  represent  a  set  of  lines  that  minimizes  that  number  of 
points  whose  distance  from  the  set  of  lines  is  greater  than  e*  yet  uses  the  least  number 
of  lines  required  to  do  so.  Hence  an  optimal  solution  of  the  partitioned  GA  does 
represent  an  optimal  function. 


4.5.3  Remarks 

1.  We  have  justified  the  claim  that  for  hm&x  >  1,  an  optimal  function  is  contained 
in  the  search  space  of  our  algorithm  and  that  our  algorithm  will  converge  to  an 
optimal  solution,  where  an  optimal  solution  represents  an  optimal  function.  We 
have  also  established  that  for  0  and  1C  large  and  critical  level  —  1  —  crit  =  1.0, 
the  solution  of  our  algorithm  will  converge  to  the  least  squares  solution  (i.e.,  the 
optimal  solution  with  h*  lines)  as  the  number  of  iterations  and  N  oo.  Hence 
our  algorithm  can  be  used  to  fit  both  robust  and  non-robust  solutions. 

2.  For  1  -  crit  <  1.0,  the  optimal  solution  to  which  the  algorithm  converges  may 
have  h'  lines  where  h'  ^  h*.  However,  a  value  of  1  —  crit  <  1.0  is  useful  when  the 
given  dataset  contains  outliers. 

3.  As  in  the  /imax  =  1  case,  for  h^ax  >  1  we  do  not  know  an  appropriate  value  for  e 
apriori.  Hence  the  proposed  algorithm  must  include  an  adaptive  procedure  which 
will  search  for  an  appropriate  value  for  e.  In  the  next  section  the  partitioned  GA 
described  above  will  be  expanded  to  include  this  adaptive  search. 

4.  A  given  dataset  of  size  N  yields  a  corresponding  estimate  oi  the  support  B  of 
the  probability  density  function  of  (X,  Y).  The  proposed  method,  by  estimating 
f{x)  via  maximization  of  FF^*^n,&,k  where  FFe*,N,Q,K  depends  upon  the  given 
dataset,  essentially  uses  B^  to  estimate  f(x).  Assuming  that  0,  /C,  and  the 
maximum  number  of  iterations  are  chosen  sufficiently  large,  given  a  sequence  of 
datsets  of  sizes  N,  N  oo,  we  may  characterize  our  estimation  of  f(x)  as  a 
maximization  over  the  corresponding  sequence  of  subspaces  Ryv  which  belong  to 
the  interior  of  B  and  approach  B  as  N  increases.  In  this  sense  our  maximization 
process  is  reminiscent  of  Grenander’s  method  of  sieves  [34]. 

The  theoretical  foundation  of  our  algorithm  has  been  established.  In  the  next  section 
we  show  the  results  of  implementing  our  algorithm  on  several  datasets  and  discuss  how 
these  results  compare  to  those  of  similar  methods. 


5  Experimental  Results 

5.1  Methods  and  Implementation 

5.1.1  Genetic  Algorithm 

We  applied  to  each  dataset  a  partitioned  GA  with  either  B  —  {2, 3, 4}  or  ?^  =  {3, 4, 5} 
(we  did  not  start  with  h  =  1  because  this  choice  for  h  was  clearly  inappropriate). 
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We  utilized  binary  coding  although  an  alternate  coding  scheme  could  have  been  used. 
Since  we  chose  values  for  0  and  K  which  were  fixed  but  very  large  (by  specifying 
la  =  8  and  Id  =  12),  following  remark  2  of  Section  4.4.3  we  implemented  a  partitioned 
GA  which  adaptively  searched  only  for  e.  We  refer  to  such  a  GA  as  a  variable  epsilon 
genetic  algorithm.  For  each  value  heH,  the  variable  epsilon  genetic  algorithm  can  be 
described  as  follows: 

1.  Set  the  global  parameters  for  population  size  M  {M  =  50),  crossover  probability 

[pc  =  0.8),  number  of  characters  for  representing  angle  la  (la  =  8),  number  of 
characters  for  representing  distance  or  offset  value  Id  (Id  =  12),  and  the  maximum 
number  of  iterations  MaxNit  {MaxNit  =  20000).  The  selection  for  M  depends 
on  the  computing  power  of  the  machine  while  la  and  Id  depend  upon  the  desired 
precision  of  our  result.  The  mutation  probability  Pm  was  varied  as  a  function 
of  the  number  of  iterations  -  see  Section  3.2  and  [29]  for  more  details.  Further 
comments  on  parameter  value  selection  can  be  found  in  Section  6. 

2.  Choose  the  critical  level  =  1  -  crit  =  percentage  of  data  points  to  fall  within  e  of 
the  final  fitted  piecewise  linear  function,  and  a  large  initial  value  for  e.  The  fitness 
value  of  each  string  (where  each  string  represents  a  function  L{0j^.^h,  kj,-,h)  € 
£“(0,/C))  is  determined  by  Equation  10. 

3.  Run  an  elitist  GA  until  either  (1)  (1  -  crit)  <  (the  maximum  fitness  value  of  the 
population) /(number  of  observations),  or  (2)  the  maximum  number  of  iterations, 
MaxNit,  is  reached.  If  (1)  occurs,  then  set  e  =  e-0.01,  Nit  =  iteration  number 
=  1,  and  restart  the  elitist  GA.  If  (2)  occurs,  report  the  function  corresponding 
to  the  string  with  the  maximum  fitness  value  as  the  final  result  for  the  given 
value  of  h. 

4.  Once  the  algorithm  has  been  executed  for  each  possible  value  of  h,  compare  the 
results  across  h  values  and  select  the  string  corresponding  to  the  overall  maximum 
fitness  value  as  the  optimum  string. 


Note  that  within  each  run  of  our  algorithm  the  value  of  h  is  fixed.  We  then  compare 
the  results  of  each  run  of  the  algorithm  for  each  value  of  h  €  ^  to  determine  the  overall 
optimal  string. 


5.1.2  Data 

The  proposed  algorithm  was  applied  to  the  four  generated  datasets  described  in  Table  1. 
Each  dataset  has  normally  distributed  noise  (denoted  Nor{mea.n,  sd))  and  was  designed 
with  specific  characteristics  -  dataset  #1:  no  outliers,  equally  spaced  knots;  dataset 
#2:  no  outliers,  unequally  spaced  knots;  dataset  #3:  outliers,  equally  spaced  knots; 
and  dataset  #4:  outliers,  unequally  spaced  knots.  The  datasets  range  in  size  from  50 
to  375  data  points. 
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Set  # 

wm 

Generating  Function 

Noise 

'  —2.0a: 

0.5a: 

,  -0.5a: +  1.0 

-1.0  <  a;  <  0.0 

0.0  <  a:  <  1.0 

1.0  <  a:  <  2.0 

Aror(0,0.25) 

H 

3 

3.88a:  +  10.44 
-1.74a: +  3.14 
,  3.77a:  -  17.25 

-3.0  <  a:  <  -1.29 
—1.29  <  X  <  3.7 

3.7  <x<  4.9 

Nor(0, 1) 

3 

60 

3 

g{x)  =  < 

'■  —0.1a: +  0.6 
-1.1a: +  1.6 
t  0.8a:  -  2.2 

0.0  <  a:  <  1.0 

1.0  <  a:  <  2.0 

2.0  <  a;  <  3.0 

Nor(0, 0.1) 

4 

50 

4 

g{x)  =  < 

(  1.3a:  +  2 

1  0.05a: +  3.125 
-0.4a:  +  3.2 
,  3.2a: -11.41 

0.0  <  a:  <  0.9 

0.9  <  a;  <  2.1 

2.1  <  a:  <  4.3 

4.3  <  a:  <  5.0 

Nor(0,0.2) 

Table  1:  Experimental  Datasets 


5.1.3  Comparison 

For  each  dataset  the  results  of  the  variable  epsilon  genetic  algorithm  were  compared 
to  the  results  of  three  other  piecewise  linear  fitting  methods: 

1.  b-spline:  a  b-spline  of  degree  2  fit  using  the  functions  bs  and  Im  in  the  software 
package  S-plus,  version  3.3,  release  1  [33].  The  knot  locations  are  equally  spaced 
by  default. 

2.  Boor  spline:  The  library  pppack  containing  functions  from  de  Boor  [13]  was 
downloaded  from  the  Netlib  repository  (http://www.netlib.org/).  The  function 
12main  was  implemented  to  fit  a  b-spline  of  degree  2  whose  knot  locations  were 
improved  by  the  subroutine  newnot.  The  number  of  iterations  was  set  at  900. 

3.  least-squares  GA:  a  genetic  algorithm  which  has  the  same  global  parameters  as 
the  variable  epsilon  GA  but  which  fits  a  piecewise  linear  function  by  minimizing 
the  sum  of  squared  error. 

All  three  methods  were  applied  to  each  dataset  for  each  value  h  E  V..  Then  for  each 
method,  the  results  for  all  h  values  were  compared  and  the  best  result  was  chosen 
as  the  method’s  overall  solution.  We  have  not  attempted  to  compare  results  under 
neural  network  based  methods  since  we  are  fitting  straight  lines  and  not  curves.  We 
did  not  compare  our  results  with  those  of  other  GA  based  methods  [2,  8,  26]  because 
the  limited  foci  of  these  works  were  incompatible  with  that  considered  here. 

All  experiments  were  run  on  a  Sun  Sparcstation  5. 

Due  to  space  constraints,  for  each  dataset  we  discuss  only  selected  results  which  demon¬ 
strate  the  relative  strengths  and  weaknesses  of  our  method.  In  the  figures  and  tables, 
the  phrase  ‘GA  Is’  refers  to  the  least  squares  GA  while  the  phrase  ‘GA  eps’  refers  to 
the  variable  epsilon  GA. 
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5.2  Selected  Results 


5.2.1  Dataset  1 

This  was  our  ‘nice’  dataset  -  no  outliers  and  a  generating  function  with  equally  spaced 
knots.  Since  there  were  no  outliers,  we  used  the  sum  of  squared  error  (SSE)  as  a 
standard  of  comparison.  The  choice  of  critical  level  for  the  variable  epsilon  GA  was 
95%.  All  four  methods  chose  the  correct  number  of  lines;  the  least  squares  GA  yielded 
the  best  fit,  followed  by  the  b-spline,  the  variable  epsilon  GA,  and  the  Boor  spline, 
respectively.  We  observe  in  Figure  5  that  only  the  Boor  spline  appears  to  be  unsatis¬ 
factory.  The  run  times  for  the  splines  were  almost  instantaneous  while  the  GAs  took 
longer  (approx.  20  minutes).  It  should  be  noted,  however,  that  with  h—Z  our  string 
length  was  60,  hence  the  number  of  possible  strings  was  1152921  *  10^^!  Despite 
this,  the  variable  epsilon  GA  yielded  an  acceptable  result  within  the  initial  5  minutes 
of  CPU  time. 


Method 

SSE 

Method 

SSE 

b-spline 

22.26 

Boor  spline 

34.23 

GAls 

17.89 

GA  eps 

23.73 

Table  2:  SSE  values  for  Dataset  1,  h=3 


Dataset  1 ,  h=3 


Dataset  1 ,  h=4 


-1.0  -0.5  0.0  0.5  1.0  1.5  2,0 


•>  X 


->x 


Figure  5:  Results  of  Dataset  1  for  h=3  and  h=4 

We  observed  that  although  the  variable  epsilon  GA  did  not  yield  the  best  fit,  it  was 
the  only  method  which  yielded  a  very  nice  fit  (SSE  =  22.80)  when  the  number  of  pieces 
was  misspecified  -  when  h  was  set  equal  to  4  although  the  generating  function  has  3 
pieces  (/i*=3)(see  Figures  5  and  6).  The  algorithm  was  almost  able  to  merge  2  of  the 
pieces  into  1.  Note  that  in  Figure  6,  where  the  result  of  applying  the  variable  epsilon 
GA  on  dataset  1  is  shown  for  h  =  4,  it  is  very  difficult  to  see  4  lines  with  the  naked 
eye.  Hence  the  variable  epsilon  GA  was  more  robust  against  misspecification  of  h  in 
relation  to  the  other  methods. 
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Dataset  1 ,  h=4 


Figure  6:  Result  of  eps  GA  on  Dataset  1  with  /i=4 


5.2.2  Dataset  2 

While  dataset  1  had  a  generating  function  with  equally  spaced  knots,  this  dataset  had 
a  generating  function  with  unequally  spaced  knots  -  the  middle  piece  was  considerably 
longer  than  either  of  the  end  pieces.  For  this  reason  the  fit  of  the  b-spline  was  quite 
poor  in  comparison  to  the  genetic  algorithms’  results.  However,  why  the  fit  of  the 
Boor  spline,  which  varies  the  knot  locations,  was  so  poor  remains  a  mystery.  The  best 
model  from  each  method  had  the  correct  number  of  lines.  The  variable  epsilon  GA 
(with  1  -  crit=95%)  yielded  the  best  model,  followed  by  the  least  squares  GA  model, 
the  b-spline,  and  the  Boor  spline. 

The  failure  of  the  least  squares  GA  to  provide  the  lowest  SSE  model  was  most  likely 
due  to  insufficient  iterations.  Again,  the  speed  of  the  spline  algorithms  was  almost 
instantaneous,  while  the  GA  algorithms  took  about  5  minutes  of  CPU  time  to  yield 
an  acceptable  result  and  about  15  minutes  of  CPU  time  to  complete  the  maximum 
number  of  iterations. 

With  this  dataset,  only  the  Boor  spline  failed  to  merge  lines  for  a  better  fit  when  the 
choice  of  h  was  large.  This  suggests  that  the  least  squares  GA  would  have  also  merged 
lines  for  dataset  1  but  perhaps  more  iterations  of  the  algorithm  were  required.  The 
b-spline  could  not  match  the  quality  of  fit  of  the  GAs  because  its  knot  locations  were 
fixed. 


Method 

SSE 

Method 

SSE 

b-spline 

610.24 

Boor  spline 

836.31 

GA  Is 

415.35 

GA  eps 

351.57 

Table  3:  SSE  values  for  Dataset  2,  h=3 
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Dataset  2,  h=3 


Dataset  2,  h=4 


-2  0  2  4  -2  0  2  4 

- >  X  - ->  X 


Figure  7:  Results  of  Dataset  2  for  =  3  and  h  =  A 

5.2.3  Dataset  3 

The  other  strengths  of  the  variable  epsilon  GA  became  evident  when  we  considered 
datasets  with  outliers.  Dataset  3  had  equally  spaced  knots  and  contained  a  cluster  of 
6  outliers.  The  outliers  made  the  SSE  criterion  useless  as  a  measure  of  fit,  so  we  based 
our  statements  about  the  quality  of  the  results  on  visual  comparisons.  Figure  8  clearly 
shows  that  the  outliers  caused  difficulties  for  all  methods  except  for  the  variable  epsilon 
GA.  The  b-spline  and  the  least  squares  GA  showed  particularly  poor  results.  However, 
the  parameter  crit  of  the  variable  epsilon  GA  made  it  robust  against  outliers.  1  -  crit 
was  set  at  90%  to  allow  for  the  outliers  -  leading,  in  this  case,  to  a  superior  fit. 

If  1  —  crit  had  been  set  at  95%,  as  was  done  previously,  the  variable  epsilon  GA  result 
would  have  been  closer  to  that  of  the  least  squares  GA  (results  omitted  due  to  limited 
space).  Even  if  the  outliers  had  not  been  clustered,  the  variable  epsilon  GA  would  have 
reached  a  superior  solution  -  see  dataset  4  for  an  example  with  non-cluster ed  outliers. 

As  with  dataset  1,  only  the  variable  epsilon  GA  yielded  a  very  nice  fit  when  the  number 
of  lines  was  misspecified  as  h  =  4  instead  of  h  =  3  (see  Figures  8  and  9).  Despite 
solution  spaces  with  sizes  of  order  ft:  1152921  *  lO’^^  for  h=S,  2®°  1208925  *  10^® 

for  h  =  4,  and  2^®“  «  1267650*  10^^  for  /i=5,  the  variable  epsilon  GA  was  able  to  yield 
a  near-optimal  result  in  about  20  minutes.  This  can  be  seen  in  Figure  9  where  the 
result  of  the  variable  epsilon  GA  for  /i  =  3  and  the  generating  lines  of  Dataset  3  have 
been  plotted  for  comparison. 


5.2.4  Dataset  4 

Our  final  dataset  combined  unequally  spaced  knots  with  the  presence  of  4  outliers 
which  were  not  clustered.  We  set  1  -  crit=92%.  From  Figure  10  it  is  clear  that 
only  the  variable  epsilon  GA  provided  a  reasonable  fit.  The  Boor  spline  and  the  least 
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Dataset  3,  h=3 


Dataset  3,  h=4 


0.0  0.5  1.0  1.5  2.0  2.5  3.0  0.0  0.5  1.0  1.5  2.0  2.5  3.0 


Figure  8:  Results  of  Dataset  3  for  /i  =  3  and  =  4 


squares  GA  were  adversely  affected  by  the  outliers,  while  the  b-spline  failed  to  capture 
the  shape  of  the  dataset.  Again,  the  robustness  of  the  variable  epsilon  GA  proved 
essential  for  a  proper  fit.  We  also  found  that  when  the  specified  number  of  lines  was 
large  (h  =  4  instead  of  h  =  3)  the  variable  epsilon  GA  was  able  to  compensate  for  this 
by  selecting  pieces  which  almost  merged. 


5.3  Conclusions 

The  above  results  demonstrate  that  for  ‘nice’  datasets  (no  outliers)  the  variable  epsilon 
GA  can  provide  a  fit  comparable  to  those  of  a  least  squares  GA  or  splines.  Although  the 
variable  epsilon  GA  did  not  always  yield  the  ‘best’  result  for  nice  datasets,  it  always 
achieved  a  satisfactory  fit  very  quickly.  When  the  number  of  pieces  was  large,  the 
algorithm  did  a  excellent  job  of  decreasing  the  effective  number  of  pieces  by  choosing 
lines  which  almost  merged.  By  increasing  the  number  of  iterations  and/or  the  string 
length,  it  is  likely  that  even  better  results  can  be  achieved. 

For  datasets  with  outliers,  the  variable  epsilon  GA  is  capable  of  yielding  results  far 
superior  to  those  of  comparable  methods.  As  with  ‘nice’  datasets,  if  too  many  pieces 
have  been  specified  the  algorithm  can  almost  merge  pieces  to  yield  a  model  with  the 
appropriate  number  of  effective  knots.  Despite  large  search  spaces  (e.g.,  on  the  order 
of  2^^  and  2®*^)  the  algorithm  can  reach  an  acceptable  solution  in  about  5  minutes  and 
a  near-optimal  solution  in  about  20  minutes.  The  variable  crit  makes  it  possible  to 
adjust  for  outliers  while  the  variables  0  and  K  provide  some  control  over  the  precision 
of  the  final  fit.  By  allowing  e  to  be  determined  adaptively,  our  algorithm  can  find  a 
result  which  (1)  satisfies  the  critical  value  and  (2)  is  closest  to  the  majority  of  the  data 
points  with  respect  to  other  solutions. 

Thus  the  effectiveness  of  the  variable  epsilon  GA  has  been  conclusively  demonstrated 
for  fitting  both  robust  and  non-robust  piecewise  linear  functions. 
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Dataset  3,  h=4 


Dataset  3,  h=3 


Figure  9:  Result  of  eps  GA  on  Dataset  3  with  /i=4;  Comparison  of  eps  GA  with 
generating  lines  on  Dataset  3  with  /i=3 

5.4  Remarks 

1.  Since  the  solution  given  by  the  variable  epsilon  GA  is  the  result  of  a  random 
process,  we  decided  to  run  an  experiment  several  times  on  dataset  5  and  examine 
the  variability  of  the  results.  Dataset  5  contained  several  outliers,  had  equally 
spaced  knots,  and  met  the  specifications  in  Table  4. 


Set  # 

N 

h 

Function 

Noise 

5 

60 

3 

g{x)  =  < 

2.5x  0.0<x  <  1.0 

-4.0a; -1-6.5  1.0  <  a;  <  2.0 

^  3.5a;  d- -8.5  2.0  <  a;  <  3.0 

Nor{0, 0.1) 

Table  4:  Dataset  5 


A  variable  epsilon  GA  was  executed  8  times  with  the  global  parameters  set  as 
in  Section  5.1.1  and  1  -  crit  =  95%.  The  results  are  shown  in  Figure  11  and 
Table  5. 

Figure  11  does  not  appear  to  contain  8  functions  because  several  results  were 
identical.  By  examining  the  ranges  of  both  the  slopes  and  intercepts  of  the  lin¬ 
ear  pieces,  the  achieved  level  of  consistency  appears  to  be  satisfactory.  Table  5 
shows  the  ranges  of  parameter  values  of  all  fitted  lines.  For  example,  the  slope 
and  intercept  of  the  first  generating  line  are  2.5  and  0,  respectively,  and  the  ex¬ 
periments  generated  corresponding  ranges  of  (2.41421,  3.02704)  and  (-0.188247, 
0.180099).  Further  research  is  needed  to  determine  a  more  precise  measure  of 
variability  of  results. 
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Dataset  4,  h=3 


Dataset  4,  h=4 


0  1  2  3  4  5 


0  1  2  3  4  5 


■>  X 


->  X 


Figure  10;  Results  of  Dataset  4  for  h  =  3  and  =  4 


Piece  #1 

Piece  #2 

Piece  #3 

Intercept 

(-0.188247,  0.180099) 

(6.64789,  6.91933) 

(-8.73460,  -7.92068) 

Slope 

(2.41421,  3.02704) 

(-4.21080,  -3.99222) 

(3.29656,  3.61354) 

Table  5:  Dataset  5  Results;  Ranges  for  Intercepts  and  Slopes 


2.  For  our  experiments  we  considered  datasets  where  h*  was  either  3  or  4  and  h  was 
between  2  and  5.  However,  the  values  for  h*  and  h  that  can  be  considered  can 
be  increased  simply  by  increasing  the  available  computational  resources.  Note 
that  when  h=A,  the  string  length  is  80  so  the  size  of  the  solution  space  is  on 
the  order  of  2®°.  Yet  we  were  able  to  achieve  near-optimal  results  in  20  minutes. 
With  greater  resources,  a  larger  solution  space  could  most  likely  be  handled  in 
approximately  the  same  CPU  time. 


6  Conclusions  and  Future  Research 

By  using  genetic  algorithms  we  have  devised  a  method  for  fitting  piecewise  linear 
functions  to  datasets  in  'R?  which  not  only  optimizes  the  number  of  pieces  but  also 
optimizes  the  knot  locations.  With  the  assumption  that  the  probability  density  func¬ 
tion  of  our  random  variables  is  symmetric,  the  above  theory  shows  that  our  method 
will  lead  to  a  piecewise  linear  function  which  ‘fits’  the  given  dataset.  However,  even  if 
we  do  not  make  this  assumption  (so  our  only  assumption  about  the  data  is  that  the 
underlying  probability  density  function  is  continuous),  our  method  will  still  yield  an 
optimal  fit  for  the  given  dataset,  optimal  with  respect  to  the  chosen  fitness  function. 

The  proposed  method  yielded  very  good  results  in  the  presence  of  noise.  The  param¬ 
eter  crit  makes  it  possible  for  our  algorithm  to  reach  a  near-optimal  result  even  in 
the  presence  of  outliers.  The  convergence  of  genetic  algorithms  has  been  shown  for 
practically  any  choice  of  parameter  values  (e.g.,  M  =  50,  Pc  =  0.8,  etc.)  although  the 
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>>  X 


Figure  11:  Example  of  Variability  of  Results 


best  choices  are  still  a  matter  of  contention.  The  formulation  of  an  optimal  stopping 
rule  is  a  subject  of  ongoing  research,  although  it  is  known  that  increasing  the  number 
of  iterations  leads  to  a  result  with  better  accuracy. 

As  mentioned  previously,  we  would  like  to  develop  a  method  for  providing  confidence 
limits  for  our  results.  It  may  also  be  possible  to  decrease  CPU  time  by  developing  a 
heuristic  for  determining  a  ‘best’  initial  value  for  e.  The  closer  e  can  be  chosen  to  its 
theoretical  optimum,  the  less  iterations  and  hence  less  CPU  time  will  be  required  by 
the  algorithm.  Given  the  speed  of  the  proposed  method  we  believe  it  is  feasible  to 
extend  the  use  of  genetic  algorithms  to  fitting  curvilinear  models  and  to  datasets  of 
dimension  greater  than  2. 
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8  Appendix 


Proof  of  Theorem  4.1: 

Let  e  >  €*,  0,  and  K,  be  given. 

Suppose  3  G  and 

{(a;,  y):y  =  kj*,i^i){x),  x  G  [zn,  ^21]}  f)  ^ 

Then  either 

{(a:,y):yG  A:j.,i,i)(a;)  +  e],  a:  G  [2:11,^21]}  fl  rect  =  $ 

or 

{(a;,y):yG  [L(^j»4,i,%*,i,i)(a;)  -  e,L(%.4,i,  rr  €  ^21]}  fl  ’’ect  =  0 

Assume  without  loss  of  generality  that 

{ix,y):ye  ^  ^  [^n,  2:21]}  f]  rect  =  0 

This  implies  that  ^  <p((^i>yi))  >  0-95,  where 


¥’((a^>y))  =  {  J 


(^>p)  ^  {(^>^)  •  y  ^  [T(%«,i,i,/jj*,i,i)(a;),L(0j.,i,i,/cj.,i,i)(a;)  +  e],  a;  G  [2:111^21]} 
otherwise 


Let  (x,y)  =  (Eili  '/^((xi,yi)))  ^  Ej(^j,yj)  where  the  sum  is  taken  over  all 
(xj,yj)  G  {(x,y)  :  y  G  [L{6j*^i^i,kj<-^i^i){x),L{6j*^i^i,kj*^i^i){x)  +e],a:  G  [^11,2:21]}. 


Define  1,1)  :  1,1,  1)  is  parallel  to  andy  =  L{6j 


*,1,1) 


Then 

{{x,y)  :  y  G  A:j.,i,i)(a;), fcj*,i,i)(a:)+e],  x  G  [2n,2;2i]}  C  ^  ^ 

Thus  G  and  since  y  =  /cj..,i,i)(x), 

{{x,y)  :  y  =  A;j..,i,i)(a:),  x  G  [2:11, 221]}  fj^ect  ^  0 


For  a  graphical  representation  see  Figure  12a.  Note  that  95%  of  the  datapoints  fall 
within  e  of  while  95%  of  the  datapoints  fall  within  e/2  of 
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Proof  of  Proposition  4.1: 

Let  LiOm, 1,U  km, I,i)  €  Pi  and  ^  >  0  be  given.  Choose  0^  :  7r/2''^  =  7r/0|  <  ^/2. 
Similarly,  choose  K.^  :  5  =  diag/{‘^^  ~  1)  =  diag/{)C^  “  1)  <  C/2-  Then  3  L{6,k)  € 
£?(0^,  )C^)  :  2.  is  satisfied.  Since  for  any  h,  h  =  1, ,  h^ax,  /C)  :  la  =  1, 2, . . . } 

and  (0,  K)  ;  Id  =  1, 2, . . . }  represent  increasing  sequences  of  nested  sets,  if  L{6,  k)  G 
£?(0^,X:^),  then  L{9,k)  G  £?(0,/C)  V  0  >  0^  and  V  Xl  >  IC^.  Hence  1  is  satisfied.  4 

Proof  of  Theorem  4.2: 

Let  ^  >  0  be  given.  Choose  0^  :  7r/2^*  =  •7r/0^  <  C/2.  Then  for  6^^^  G  •  •  •  >  tt};  < 
%+i)  ^  b  i  =  1, . . .  ,  0$  -  1},  we  have  1 0  -  1<  (/2, 

I  %)  -  %(2)  l<  C/2,  -  -  •  ,  I  -  TT  |<  ^/2.  So  given  any  L{9m,i,u  km,i,i)  e  Pi  we  can 
choose  0^  so  that  3  G  {^, . . .  ,7r}  :  1  9m,i,i  -  |<  C/2- 

For  any  angle  G  [^,  n  =  1, . . .  ,  0^  -  1  ,  the  corresponding 

P?(n)  e  <  %nO  ^  T'«(n,)  <  ^^“5-  Find 

1/  =  max  {sup  n  =  1, . . .  ,  0^  -  1}} 

kn) 

Choose  JC^  :  z^//C^  <  C/4  and  fC^  =  2^  for  some  yeU.  Then  given  any  L{9m, i,i,  km,i,i)  € 
Pi  we  can  choose  so  that  3  G  (0, . . .  ,  2^«  —  1}  :  |  —  Pm,i,i  \<  C/2- 

Hence  given  any  C  >  0  L{9m, i,i,  km, i,i)  €  Pi  we  can  find  0^  and  /C^  so  that 

3  L(%),  e  £?(0e,/C^)  :|  9  -  9m,i,i  \<  C/2  and  |  p^,i,i  -  p  |<  C/2. 

If  L{9^,,^,k^J  G  £?(0e,X:^),  then  G  £5(0, /C)  V  0  >  0^  and  V  AC  >  AC^. 

Hence  1.  and  2.  are  satisfied.dfk 

Proof  of  Proposition  4.2: 

Define  Si  =  U£i  £i(0i,/Ci).  Note  Pi  C  5i.  Note  that  <Si  C  7{  so  Pi  C  Since  we 
are  requiring  our  proper  lines  which  pass  through  rect,  the  proper  lines  in  Pi  =  the 
proper  lines  in  Ti- 4 

Proof  of  Theorem  4.3 

Note  that  Ci{Qi,)Ci)  C  £5(0,', Xli')  W  i'  >  i.  Hence  Ad  Q  Aev  y  i'  >  i.  Thus 
.Aeiim  =  limi_^.ooAi  exists  [30].  But  AeHm  is  the  set  of  optimal  lines  in  Si  where 
Pi  ^  «5i  C  7i.  By  Proposition  4.2,  Aiim  =  «4e- 

Proof  of  Theorem  4.4 

Suppose  3  £(^j*,i,i,  A:j.,i,i),  £(0j*«,i,i,  A:p*,i,i)  e  .4ej^(1.0)  : 

£(%*,i,i,  ^i*,i,i)  ^  £(%**, 1,1,  ^j**,i,i)-  Then  for  all  possible  values  of  (x,y) 

(x,y)  G  {{x,y)  :  y  G  [L(0j.,i,i,  A:j*.i,i)(a:)  -  e{  o,  L(6’p,i,i,  kj*^i^i){x)  +  q]} 

and 
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(x,y)  G  {{x,y)  :  y  G  -e^  o,  (,]}  . 

Recall  that  for  h*  =  1,  the  support  of  a\{x,y)  was  defined  as  B  where 

B  =  {(:£,  v):y€  -  £o,  +  «ol}  for  x  €  lz„,  Zji] 

As  iV  -)■  oo,  Cl  0  Co  since  e*  g  is  minimal  and  is  continuous  for  all  e.  Hence 

{{x,y)  :  y  =  {{x,y)  :  a;cos6>i  +  ysin^i  =  dj  and  {{x,y)  : 

y  =  ^  {{x,y)  :  a; cos + 1/ sin 6>i  =  dj  Thus  = 

See  Figure  12b  for  a  graphical  representation  of  Theorem  4.4. 


(a)  (b) 

Figure  12:  Graphical  representations  of  (a)  Theorem  4.2  and  (b)  Theorem  4.4 

Proof  of  Theorem  4.5: 

Given  the  continuity  of  FFe,N  we  know  that  as  iV  ^  oo,  for  each  N  large  we  may  find 
^iN)  £2W)  csiv,  0  <  ciiv  <  e2N  <  C3W)  such  that 

•  AFFe^j^ kj*,i,l)eiN,N)  <  0.95 

•  A:j*,i,i)e2jv,Jv)  >0.95 

•  AFFe3jy,yv'(T(dj»,i,i,  >  0.97 

Note  that  at  0  £(eiN)  while  is  not  minimal,  i.e.,  there  exists  some 

e  <  esjv  :  ^  0-  Hence  for  our  stated  goal  (see  Eqn.  11)  we  are  interested  only  in 

e2N-  For  each  N  there  may  be  infinitely  many  such  e2N-  For  a  given  N  and  e  let  one 
such  e2N  be  Then  we  conclude 

liminfAr^oo^FFe.^,Ar(L(6»j.,i,i,A:j.,i,i)e*^,Ar)  >  0.95 

for  appropriate  0  and  JC.  4)k 
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